PSR  Report  2628 


AD“A286  887 


» 

NoiseProp:  A  Dynamic  60  kHz  to  30  MHz  Atmospheric 
Radio  Noise  Model 

ft 


Ronald  M.  Bloom 

David  M.  Crandall  • 

Chris  R.  Warber 


March  1996  »  • 

it 

Ajpprovec  Im  fit.'j'  .ei&ojMS 

Final  Report 

Contract  N00014-94-C-0053 


Prepared  for 

Office  of  Naval  Research 
Ballston  Tower  One 
800  North  Quincy  Street 
Arlington,  VA  22217-5660 


9  6  6  2 1  5  1  7 


96-01314 


ft 


^  ISACIFIC-SiERRA  RESEARCH  CC?RRORATICBI\l 

~  ~  2r01  28th  Street*  Santa  Monica.  CA  90405  •  (310)  314-2300,  FAX  (310)  314-2323 

I 

F-wtuL'crsD  1 


•  •  • 


•  «  # 


e 


» 


27.  Abstract: 


The  combined  progation  and  lightning  source  algorithm  are  contained  in  a  \A/indow3-95  based  graphical  user  interface 
she'!  called  MoiseProp.  This  program  uses  as  data  either  the  original  satellite-based  GLO  maps  or  the  dynamic, 
meteorologically  based  DGLO  maps  as  they  become  available. 


28.  Classification  of  Abstract: 

Unclassified 


29.  Limitation  of  Abstract; 

Unlimited 


31.  Required  Peripherals: 

32.  U  of  Physical  Records: 

33.  #  of  Logical  Records: 

34.  U  of  Tracks: 

35.  Record  Type: 

38.  Color: 

37.  Recording  System: 

38.  Recording  Density: 

39.  Parity: 

40.  Playtime: 

4’i.  Playback  Speed;  42.  Video  43.  Text;  44.  S 

46.  Other:  NOISPROP  requires  WINDOWS95  and  0  MB  of  hard  disk  space 


43.  Text;  1 44.  Still  Photos:  1 45.  Audio: 


47.  Documentation/Suppiementai  infonnation: 

Also  includes  hardcopy  report:  NoiseProp:  A  Dynamic  60  Khz  to  30  Mhz  Atmospheric  Radio  Noise  Model 


1 43.  Point  of  contact/Teiephone  #; 

I  Office  of  Naval  Research;  703-696-6993 


DTIC  FORM  530  JUN  94 


Security  Classification  of  this  Pago 


:-lAi 


PREFACE 


i*) 


This  report  describes  a  broadband  (60  kHz  to  30  MHz)  electromagnetic  background  noise  model 
(NoiseProp),  The  mcdel  is  conceptually  similar  to  the  PSR  longwave  noise  model  LNP,  re¬ 
leased  in  1991.  Global  distributions  of  lightning  activity  divided  into  .seasonal  and  diurnal  maps 
are  used  to  determine  a  set  of  elemental  noise  transmitters.  Ilie  power  radiated  by  each 
transmitter  is  proportional  to  the  lightning  flash  rate  at  that  location  and  to  an  empirically  de¬ 
termined  energy  spectrum  (which  varies  roughly  as  f'^).  The  radiated  power  is  then  propagated 
by  a  variety  of  propagation  algorithms  (LF,  MF,  HF)  to  the  receiver.  The  overall  noise  power  at 
the  receiver  is  taken  to  be  an  incoherent  sum  of  power  propagated  from  each  significant  noise 
transmitter. 

The  report  additionally  describes  methods  developed  which  will  allow  certain  up-to-date  and 
forecasted  weather  data  to  be  converted  to  lightning  activity  (flash  rat;)  maps.  In  this  way,  the 
NoiseProp  model  can  be  made  to  predict  amiospheric  noise  dynamically  rather  than  simply  long 
term  median  values. 

A  Windows-95  based  graphical  user  interface  for  the  NoiseProp  has  been  developed  and  is  de¬ 
scribed  in  the  Appendix. 
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1.0  INTRODUCTION 
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In  1991  Pacific-SieiT!!  Research  (PSR)  Corporation  completed  a  PC-basea  lonj;  wave  noise 
prediction  model  iLNP)  [Warber  and  Field,  1991  ].  LNP  calculates  worldwide  median  narrow 
band  noise  power  in  the  ELF  (50  to  300  Hz)  and  VLF  (10  to  60  kHz)  bands  by  summing 
radiated  power  from  worldwide  distributions  of  lightning  flashes.  Tlie  lightning  flash  data  are 
contained  in  seasonal  and  diunial  global  lighting  occurrence  (GLO)  maps  synthesized  from  two 
years  of  the  Japanese  L'JS-b  satellite  (1978  to  1980)  HF  sensor  lightning  observations.  LNP  was 
designed  as  a  physics  based  alternative  to  the  standard  CCIR  empirical  global  noise  radio  model 
[CCIR,  1988],  comparable  to  but  more  sophisticated  than  the  WGL  noise  model  [Maxwell  ct  al., 
1970]. 

The  work  described  in  this  report  consists  of  extensions  to  the  predictive  methodology  of  LNP 
applied  to  the  joint  problem  of  1)  LF,  MF,  HF  radio  noise  prediction  (60  kHz  to  30  MKz)  and  2) 
development  of  a  capability  to  supplement  the  original  long-term  historical  GLO  maps  with 
short-term  dynamic  global  occurrence  data  (DGLO)  based  on  certain  relevant  meteorological 
data. 

•  Work  detailing  the  correlation  between  lightning  flash  rate  and  certain  commonly  available 
meteorological  parameters  (and  plans  for  the  association  of  such  data  into  the  LNP  noise 
prediction  system)  was  done  under  DNA  sponsorship  and  is  described  completely  in  PSR 
Report  2605  [Warber  and  Sinclair,  1995].  By  augmenting  the  physical  basis  of  the  noise 
model  in  this  manner  it  will  be  possible  to  predict  radio  noise  for  short  or  near  real  time 
meteorological  conditions,  I'he  existing,  widely  u.sed  models  (e.g.,  CCIR)  rely  strictly  on 
long-term  averages  and  cannot  adjust  to  specific  weather  (hence,  lightning  Hush  rate) 
conditions. 

•  The  existing  LNP  algorithms  propagated  ELFWLF  RF  energy  in  an  earth-ionospheric 
waveguide  quasi-exponential  mode  sum  along  each  great  circle  path  from  lightning  source 
region  to  receiver.  {Actually  all  lightning  activity  in  a  10°  X  10°  ceil  is  referenced  to  an 
equivalent  noi.se  tran.smitter  (ENT)  at  the  center  of  ma.s.s  of  the  tTll.j  Tlie  waveguide  mode 
formulation  |  Budden,  1966]  is  impractical  for  frequencies  above  ~100  kHz,  as  large  numbers 
of  discreie  inodes  must  be  summed.  The  .search  for  discrete  modes  at  higher  frequencies  also 
is  hampered  by  the  tendency  of  the  modal  parameters  (propagation  constants)  to  coalesce  as 
the  wave  frequency  is  raised.  Our  propagation  algorithms  will  therefore  be  different  from 
LNPs.  They  are  a  combination  of  ray-methods  with  some  diffractive  corrections  and  purely 
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empirical  methods  used  in  the  transitional  upper  MF  lower  HF  band  where  the  closeness  of 
the  gyro-frequency  is  a  crucial  determinant  of  propagation. 

•  Additionally,  the  original  lightning  stroke  spectra  used  in  LNP  have  been  extended  above  60 
IcHz  using  the  available  empirical  data. 

2.0  NOISEPROP  PROGRAM 

The  combined  propagation  and  lightning  source  algorithm  are  contained  In  a  Windows-95  based 
graphical  user  interface  shell  called  NoiseProp.  The  NoiseProp  program  uses  as  data  either  the 
original  satellite-based  GLO  maps  or  the  dynamic,  meteorologically  based  DGLO  maps  as  they 
become  available.  NoiseProp  allows  the  user  to  select  one  or  .several  locations  using  an 
interactive  map.  Subsequently,  radio  noise  in  a  given  bandwidth  as  a  function  of  local  time  at  a 
fixed  frequency  (or  ate  .several  frequencies  at  a  fixed  time)  i.s  computed  and  displayed.  The 
noise  level  is  normally  the  median  (50^  percentile).  Noise  exceeded  an  arbitrary  (0.0001  to 
99.999)  percentage  of  time  can  also  be  computed.  The  noise  is  in  units  of  dB  relative  to  kTb 
where  k  is  Boltzmann's  constant,  T  is  an  absolute  reference  temperature,  b  is  bandwidth  in  Hz 
measured  in  a  short  vertical  antenna  over  a  conducting  ground  plane.  This  is  identical  to  the 
CCIR  322-3  [  CCIR,  1988]  scaling  and  was  used  in  the  final  calibration  of  the  NoiseProp 
program. 

3.0  PROPAGATION  OP  RF  NOISE  FROM  GLOBAL  DISTRIBUTION  OF  LIGHTNING 
ACTIVITY:  60  kHz  TO  20  MHz 

As  in  the  LNP  ELF.  VLF  prediction  code  Noi.seProp  divides  the  globe  into  10°xl0°  sectors. 

Each  10“  X  10°  sector  is  regtuded  as  an  equivalent  noi.se  transmitter  (ENT).  The  long  term 
.median  lightning  flash  rates  computed  from  the  ISS-b  satellite  GLO  maps  are  sorted  by  l(l“xlO° 
sectors.  Tlie  accumulation  of  individual  noise  transmitters  is  combined  into  a  single  effective 
radiator  —  the  ENT— and  positioned  at  tiie  center  of  mass  of  the  sector.  Typical  flash  rate  data  is 
shr)wn  in  Figure  3- 1 .  The  power  contributions  from  all  worldwide  ENT;;  are  then  summed  to 
yield  the  total  median  noi.se  power  (at  a  particular  receiver  location).  This  procedure  is  the  .same 
whether  GLO  data  or  short  term  DGLO  data  are  used.  If  DGLO  data  updates  are  only  available 
over  limited  geographical  areas,  both  GLO  and  DGLO  maps  will  be  used  with  the  original  long 
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Figure  3-1.  Typical  global  lightning  flash  rate  distribution.  From  Jiine-July-August,  2200-0200  UT  ISS-b 
satellite  data. 


term  GLO  maps  rc’placing  the  DO)  ,0  maps  where  the  latter  are  nonexistent.  This  method  is 
illustrated  in  Figure  3-2  and  the  summation  can  be  written 

N(x,y)  =  ^P,(x,y)p,n,  ,  (3.1) 

I..I 

where  the  summation  is  over  all  ENTs  and 

N  N  is  the  ttrtal  noise  power  at  receiver  locution  (x,y)  in  a  1  Hz  bandwidth 

Pj  is  the  power  in  watts  of  the  i*'*  ENT 

rij  is  the  number  of  flaslies  (in  1  s)  represented  by  the  i‘h  tiNT 

Pj(x,y)  is  the  propagation  factor  (loss)  for  the  path  from  i*''  ENT  to  the  location  (x,y). 

The  piupdgation  factor  Pj(x,y),  in  principle,  depends  on  the  radiowave  frequency,  the  great  circle 
distance  between  ihe  ih>  ENT  and  the  receiver  position  (x,y),  and  the  complete  geophysical 
description  of  the  propagation  path;  namely,  the  ionospheric  conductivity  height  possible,  the 
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F-igure  3-2,  Noiso  model  methodology:  Equivalent  Noise  Transmitters  (ENTs)  and 
propagation  paths  to  a  receiver. 

ground  conductivity,  and  the  geomagnetic  field  elements  specified  at  each  position  along  the 
patli  (at  a  particular  time  and  season).  P,  repre.sents  the  propagation  of  one  watt  from  the  i^^’ 
ENT  and  tP;  .pj]  repre.sents  the  propagation  of  the  energy  in  one  flash  Ifom  the  i*''  ENT. 

In  the  ELF,  VLF  regime  (wavelengths  ~  10  to  1000  km)  the  path  is  subdivided  into  segments 
whose  ionospheric,  geomagnetic,  and  ground  properties  can  be  considered  constant  (or  insigni¬ 
ficantly  varying  on  the  scale  of  a  wavelength).  Each  such  segment  then  provides  boundary 
conditions  for  the  solution  of  an  eigenvalue  problem  for  the  propagation  constant  (mode 
constant)  in  the  model  expansion  for  the  transverse  electric  field  in  the  earth-iono.sphere 
waveguide  |  Davies.  1990;  Budden,  1961;  Warber  and  Field,  1991 1.  Modal  expansions  on  each 
patli  segment  are  matched  so  as  to  satisfy  certain  boundary  conditions  (mode  coupling)  and  the 
propagation  is  then  complete.  At  ELF  and  VLF  only  a  few  inodes  lue  needed  for  an  adequate 
expansion  of  the  field.  Mode  constants  for  various  carth/ionospheric/geomagnetic  parameter 
values  are  precomputed  and  interpolated  as  needed,  making  the  propagation  algorithm  fairly 
rapid  and  efficient. 

ITie  relative  simplicity  of  the  modal  expansion  breaks  down  above  VLF  and  we  must  use 
alternate  approaches  to  the  solution  of  the  propagation  factor.  At  the  MF  and  above  (~100  kHz 
to  30  MHz)  the  wavelength  is  such  that  the  ionosphere  can  be  said  to  be  a  slowly  varying 
medium  witli  respeci  to  a  wavelength.  T'he  methods  of  geometrical  optics  arc  applicable  for  the 
determination  of  signal  power  by  multiple  ionospheric  and  ground  hops  (see  illustration  in 


•  « 


9 


Figure  3-3).  Geometrical  optic.s  and  ray  tracing  vnethod.s  for  HF  propagation  have  been 
developed  for  several  decades  and  exist  in  various  degree.s  of  sophistication.  They  range  in 
complexity  from  simplified  engineering  formulations  originally  used  with  charts  and 
nomographs  |CRPL,  19481  to  full-blown  shooting  methods  which  trace  rays  through  ionospheric 
profiles  with  realistic  profiles  and  horizontal  gradients  (Hatfield  and  Smith,  19871.  Propagation 
losses  are  eitlier  accumulated  using  .seini-empirical  rules  at  each  ionospheric  hop  (in  the  simpler 
methods)  or  accrued  by  analytically  or  rumerically  integrating  the  ionospheric  refractive  index 
along  the  ray  path  as  it  is  traced  (in  .'he  sophisticated  methods).  Naturally,  what  is  gained  in 
accuracy  with  the  soph.bsticated  ray  tracing  methods  is  usually  lost  in  efficiency.  These  methods 
must  perfonn  time  consuming  iterations  to  find  the  exact  ray  solutions  connecting  a  tran.smitter 
and  receiver.  Furthermore,  the  accuracy  in  such  a  method  may  be  greater  than  the  accuracy 
intrinsic  to  the  underlying  ionospheric  model. 
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Figure  3-3.  Wave  guide  propagation  for  ELF-VLF  signals:  ionospheric  and  geomagnetic  properties  vary 
significantly  within  one  wavelength.  Geometrical  optics  propagation  for  LF-MF-HF  propaga¬ 
tion:  ionosphere  can  be  considered  slowly  vaiying  on  the  scale  of  one  wavelength. 
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A  commonly  used  approximate  method  for  communications  purposes  is  called  the  virtual  ray 
trace  method.  In  this  method,  the  propagation  path  is  divided,  a  priori,  into  a  certain  number  of 
hops  (depending  on  the  path  length  and  the  iono.spheric  profile  at  one  or  more  positions  along 
the  path).  On  each  segment  an  appropriate  takeoff  angle  and  equivalent  reflection  height  is 
determined  by  an  empirical  or  analytic  relationship  between  the  hop  length  and  ionospheric 
profile  at  the  hop  midpoint.  Assumption  of  a  special  (e.g.,  parabolic)  electron  density  profile 
facilitates  analytic  calculation  [Barghausen  et  al.,  1969|.  The  resulting  set  of  independent  hops 
are  then  combined  into  one  geomec  ically  self-consistent  mode  by  a  puss-range  adjustment 
algorithm  which  guarantees  equal  takeoff  angles  on  all  hops. 

After  extensive  investigation  of  the  existing  alguritlims  for  HF  (2  to  30  MHz)  propagation,  we 
developed  a  hybrid  method  which  is  a  virtual  ray  trace  method  based  on  standard  ionospheric 
maps,  but  is  not  burdered  by  the  communications  oriented  overhead  of  the  canonical  programs, 
e.g.,  lONCAP.  Indeed,  ia  NoiseProp  application,  dozens  or  even  hundreds  of  noise  transmitters 
(ENTs)  are  propagated  over  a  random  assortment  of  paths  and  we  would  expect  tliat  subtle 
differences  between  several  different  propagation  methods  be  washed  out  in  tlie  overall 
summation  of  energy.  Therefore,  computationai  expediency  can  be  regarded  as  a  priority  of  the 
first  order.  (Even  as  such,  the  method  ultimately  chosen  is  quite  computationally  iutencive  over 
a  grid  of  many  receivers.) 

The  method  used  in  the  LF-MF  regime  (60  to  900  kHz)  is  also  a  skywave  wave  hop  method.  It 
differs  from  that  used  at  HF  chiefly  in  a  different  ionospheric  los;:  rrodel  (empirical  reflection 
coefficients)  and  diffractive  corrections  required  at  the  longer  LF-MF  wavelengths  (but  not  at 
HF).  Between  90C  kHz  and  2000  kHz  a  third  empirical  propagation  algorithm  was  chosen.  In 
this  frequency  range,  the  magnitude  and  direction  of  the  earth's  magnetic  field  becomes  a  more 
.sensitive  determinant  of  propagation  then  elsewhere.  Accordingly  a  propagation  algorithm 
specific  to  this  upper  MF  band  is  employed  in  NoiseProp. 


■■K) 
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4.0  EXTENDING  LIGHTNING  SOURCE  SPECTRA  ABOVE  60  kHz 


The  original  LNP  source  model  consists  of  a  hypothetical  ENT  located  at  the  center  of  each 
i0°xlU‘^  .sector  of  the  earth.  The  power  radiated  in  a  1  Hz  bandwidth  by  an  ENT  is 


P  =  ESD  «  FR  (watts)  . 


(4.1) 


The  ESD  is  Ae  generic  flash  energy  spet-tral  density  in  units  of  Joules/Hz  (or  equivalently 
watts/HZ“)  and  FR  is  the  flash  rate  for  tne  particular  10°xl()°  cell  in  units  of  sec*.  The  flash 
rate  is  determined  directly  from  the  long  term  satellite  based  GLO  maps  (or  the  short  teira 
meteorology  based  dynamic  DGLO  maps  us  they  become  available).  The  ESD  must  be  derived 
from  various  lightning  spectra  measurements  reported  in  the  literature. 

Typical  lightning  spectra  measurements  are  made  at  distances  of  10  to  50  km  from  a  flash  and 
are  detected  either  by  wide  band  receiver.s  or  banks  of  several  narrowband  receivers  (tuned  to 
different  frequencies)  [Uman,  1987;  Weidmanand  Krider,  1986;  Nanevicz,  1987;  Willet  et  ah, 
1990].  When  banks  of  narrow  band  receivers  are  used  the  spectrum  of  a  flash  can  be  plotted 
directly,  but  the  resulting  spectrum  will  be  a  sum  of  all  flash  constituents  spectra  (strokes)  which 
occur  within  the  averaging  time  of  the  naiTow  band  filtering;  stroke  spectra  within  a  flash  cannot 
be  individuated  by  narrow  band  measurements.  Wide  band  receivers  coup'cd  to  digital  fourier 
analysis  systems  can,  on  the  other  hand,  be  used  to  find  separately  the  spectra  of  the  flash 
subprocess  or  strokes;  namely,  the  return  stroke  (RS),  stepped  leader  (SL),  in-cloud  (IC),  etc. 
[Warber  and  Field,  1991]. 

Mathematically, 

S(t)  =  Xsi(t-'C,)  (V/m),  (4,2) 

t-l 

where  S(t)  is  the  measured  electric  field  as  a  function  of  time  of  the  entire  flash  and  the  Sj  are  the 
constituent  strokes  within  the  flash.  The  Xj  are  the  times  at  which  the  i*^  stroke  commences.  We 
assume  that  the  strokes  within  a  flash  do  not  overlap  and  that  the  Xj  are  otherwise  randomly 
distributed.  It  then  follows  that  the  flash  spectrum  is  given  by 

l|S(f'f  =X|s,(f)||  (V/m)’ /Hz'  (4.3) 

(=1 

where,  by  definition 

S(f)^  js(t)e'’"''dt  ,  (4.4) 

Narrowband  measurements  tend  to  reproduce  values  of  ||S(f)||'  directly,  whereas  a  posteriori 
analysis  of  digitized  wideband  signals  allows  separate  reconstruction  of  the  individual  stroke 


k} 


spectra  ||S(f)f  i=l, N.  Corresponding  to  each  sj  is  an  energy  spectrum  (ESDj)  and  the  ESD 
lor  the  entire  flash  is  the  sum  of  the  ESDs  of  the  constituent  strokes. 


We  show  next  that  at  medium  and  high  frequencies  (>100  kHz)  the  energy  spectral  density  ESD 
is  related  to  the  flash  spectrum  by 


ESD  = 


30 


watts  /  Hz  ^ 


(4.5) 


where  d  is  the  distance  from  flash  to  the  observation  point. 

In  a  narrowband  receiver  the  i '  ponse  due  to  radio  frequency  energy  from  a  lightning  flash  is 
seen  as  an  oscillation  at  the  bandpass  center  frequency  f  which  is  amplitude  modulated  by  a 
slowly  varying  envelope  e(t). 

E(t)  =  e(t)»sin(2mft  f  (t))dt  V  /  m  (4.6) 


-«) 


Now  letting 


S(f)  =  J  E(t)e ’"*^'dt 


(4.7) 


us  above,  we  have 


5(f)||-df  =  jjF/(t)dt 


(4.8) 


j  E^(t)dt  =  je"(t)sin^(27tft  +  ct))dt  =  — Je’(t)dt-—  j'e'(t)cos(47cft +  (j))dt  (4.9) 


SO 


jE-’{t)dt  = 


(4.10) 


.since  the  last  integral  is  vanishingly  small  when  e(t;  is  slowly  varying.  Hence, 


S 


•  •  • 


•  • 


m 


m 


a 


(4.11) 


J||S(r)||'cif  =-^-Je-’(t)dt  . 


Letting  BW  denote  the  receiver  bandwidth: 


||S(f)ir  = — [e'(t)dt 
'  "  4 .  BW  ■' 


(4.12) 


Now,  at  MF  and  HF  it  is  appropriate  to  consider  the  lightning  flash  as  an  isotropic,  randomly 
oriented  radiator  at  approximately  the  cloud  base  height  (2  to  7  kxn)  [Uman,  1987;  Horner,  1964; 
Kotaki,  1984].  If  e(t)  is  the  instantaneous  envelope  in  V/m  measured  in  a  short  vertical 
monopole  at  ground,  the  instantaneous  power  flux  will  be 


e‘ 


2Z„  240jt 


-W/nri 


(4.13) 


We  then  integrate  the  isotropically  distributed  flux  over  a  hemisphere  of  radius  d  (d  is  the 
horizontal  distance  from  source  to  observer  (d  ~  10  to  50  krn)  to  obtain  total  instantane  ms 
source  power 

e'd' 


P  = 


120 


(4.14) 


The  total  energy  in  the  flash  (in  the  measurement  passband)  is  E  =  J  pdt  Hence 

1 2  I  2 

E= - |c^(t)dt  =  —  •  BW»|S(f)|p  watt -sec  ,  (4.15) 


u.sing  Eq.  (4.12)  above. 

Consequently,  the  mican  power  (using  Eq.  xx  above)  over  the  time  Tfsec)  of  the  flash  duration  is 
(in  bandwidth  BW); 


(P)  = 


d- 


nw||S(f)|p 

3()T 


watts 


(4.16) 


and  the  energy  spectral  density  (ESD) 


if) 


ESD- 


d^  .||S(ff 


30 


watts  / 


v«) 

(4.17) 


In  the  LNP  formalism,  below  60  kHz,  the  ESD  is  built  up  from  a  weighted  sum  of  dominant 
stroke  spectra||Si  (f  )||^ .  These  spectra  Sj(f)  are  found  for  return  .strokes  and  certain  in-cloud 

radiators  in  articles  by  Weidman  and  Kxider  [1986];  Willet  et  al.  [1990],  Uman  [1987];  and 
Greifinger  [1989].  We  have  extended  this  decomposition  of  the  overall  flash  ESD  from  60  kHz 
up  to  1  MHz,  utilizing  the  same  weighting  as  in  LNP.  In  LNP,  a  typical  flash  will  consist  of 
four  return  strokes.  20  vertical  in-cloud  strokes,  and  20  horizontal  in-doud  strokes  [Warber, 
1991],  Tlie  number  of  strokes  per  flash  is  an  adjustable  parameter  of  the  model.  Above  about  1 
MHZ,  however,  the  simple  description  of  a  flash  as  consisting  of  a  small  number  of  recurring 
basic  stroke  types  becomes  inadequate.  Above  1  MHZ,  a  narrowband  receiver's  response  to  a 
lightning  flash  is  significantly  enhanced  by  a  composite  of  hundreds  or  tJiousands  of  rapid 
intercloud  discharges  which  are  not  strong  radiators  at  VLF  [Nanevicz,  1987;  LeVinc,  1987], 

Hence,  at  HE,  the  combined  spectra  of  return  strokes,  leader  pulses,  and  tlie  few  other  specific 
stroke  types  that  have  been  separately  analyzed  in  the  literature  cannot  reproduce  the  spectral 
amplitude  and  slopes  seen  in  typical  n;uTowband  receiver  measurements.  We  have  therefore 
used  an  empirical  amplitude  spectrum  for  the  entire  flash  at  frequencies  above  1  MHZ  [Homer, 
1964;  Kotaki,  1983,  1984).  Figure  4. 1  shows  the  extension  of  the  return  stroke  and  bipolar  (in¬ 
cloud)  stroke  spectra.  The  curves  below  100  kHz  are  identical  with  those  used  in  the  original 
LNP  model.  Between  100  kHz  and  MHz,  tnc  return  stroke  spectrum  and  the  bipolar  spectrums 
are  interpolated  so  as  to  Join  smoothly  with  the  return  stroke  .spectrum  shown  in  Willet  et  al., 
[1990]  above  1  MHz,  The  specti'a  in  Weidman  and  Krider  [1981,  1986]  indicate  that  the  bipolar 
in-cloud  stroke  spectra  are  essentially  coincident  witii  the  return  .stroke  at  higher  frequencies. 
Above  I  MHz  the  estimated  flash  spectrum  of  Horner  ]  1964]  shown  by  tlie  straight  line  in 
Figure  4-1  is  used  instead  of  the  LNP  weighted  combination  of  sti’oke  spectra. 


Due  to  the  spar  sity  of  published  complete  flash  characterizations  above  V  LF-LF,  the  spectral 
models  in  the  NoiseProp  model  are  inevitably  a  source  of  uncertainty  in  the  final  results,  namely 
the  noi.se  at  a  given  receiver  location.  A  posteriori  calibrations  of  the  model  output  as  a  function 
of  frequency  must  be  made  (see  Section  1 1,0  under  calibration;. 


9 


9 


9  • 


O 


10  10^  103  10**  10® 

I  Freq.  (kHz) 


Warber  data  Interpolation  Willet  data 

Figure  4-1 ,  Extended  ENT  source  spectra.  Dimensions  are  that  of  Energy  Spectral  Density  (ESD)  in 
Watts/Hz/Hz  (Joulea/Hz).  Return  stroke  spectrum  is  interpolated  between  the  LNP 
spectrum  and  the  spectrum  of  Willet  [1990],  The  vertical  bi-polar  spectrum  is  smoothly 
interpolated  so  as  to  join  that  of  the  Willet  [1990]  Return  Stroke.  The  Horner  [1 964]  com¬ 
plete  narrow-band  flash  spectrum  is  also  shown. 


5.0  SUMMARY  OF  LF,  MF  PROPAGATION  METHODS 


Low  frequency  radio  waves  (30  to  300  kHz)  are  reflected  at  the  lower  side  of  the  iono.spheric  D- 
layer  (-TO  kni,  day;  -90  km,  night)  and  may  propagate  strongly  over  distances  up  to  1000  to 
5000  km.  Above  60  kHz  (the  upper  limit  of  the  LNP  propagation  capability)  it  becomes 
convenient  to  represent  propagation  as  occurring  along  geometric  hops  defined  by  a  set  of 
ionospheric  and  ground  reflections  jJohler,  1962;  CCTR,  1990].  This  is  referred  to  as  the 
wavehop  method. 

Medium  frequency  (300  to  3000  kHz)  propagation  by  U  and  E-layer  reflection  is  strongly 
affected  by  the  geometry  of  the  eartii's  magnetic  field  since  the  electron  gyro  frequency  (~0.8  to 
1 .6  MHz)  always  falls  somewhere  within  this  band.  Moreover  D-layer  absorption  of  MF  signals 
is  quite  pronounced.  During  daylight,  distance  skywave  MF  signals  are  \  ery  highly  attenuated 


and  a  local  ground  wave  signal  will  dominate.  Nighttime  D-layer  absorption  is  .small,  allowing 
long  distance  propagation.  The  complexity  of'  a  purely  physical  propagation  model  at  MF  is 
such  that  empirical  or  semi-empirical  models  are  u.sually  employed  [Davies,  1990;  PoKempner, 
1980;  CCIR,  1992bl. 

6.0  NOISEPROP  LF  PROPAGATION  METHOD 

Following  CCIR  265-7  [CCIR,  1990)  and  Johlerll96]  I  we  have  developed  a  wavehop  method 
that  allow.s  calculation  of  tire  radiated  field  from  a  lightning  source  at  frequencies  from  60  kHz 
up  to  the  middle  MF  band  (500  to  900  kHz).  The  wavehop  method  is  a  geometrical-optics 
approximation  with  added  diffraction  corrections  to  account  for  the  sphericity  of  the  earth  and 
ionosphei'e  [Wait,  I960].  The  algorithm  divides  the  path  between  the  ENT  source  and  the 
receiver  into  a  sequence  of  hops.  Tliis  is  accomplished  by  an  exact  geometric  solution 
connecting  an  elevated  source  and  a  ground  based  receiver  by  up-going  and  down-going  ray 
paths,  allowing  for  tire  crossing  of  the  day-night  terminator.  Day  and  night  reflections  are 
assumicd  to  occur,  to  a  first  approximation,  at  70  and  90  km,  respect'  /ely,  at  lower  frequencies 
(D-layer  reflection)  and  at  105  and  1 10  km  at  higher  frequencies  (E-layer  reflection).  Figure  6-1 
shows  a  schematic  of  a  multihop  day  to  night  path.  Tliese  reflection  heights  are  interpolated  as 
functions  of  wave  frequency  between  60  and  900  kHz. 

Associated  with  each  ionospheric  reflection  point  is  a  2  X  2  matrix  of  reflection  (or  transmission) 
coefficients  describing  the  ratios  of  reflected  anil  incident  w  aves  below  the  ionosphere  of  two 
ortliogunal  polarizations.  The  magnitudes  of  the  reflection  coefficients  describe  the  attenuation 
of  field  strength  incurred  in  an  ionospheric  reflection.  Ilic  reflection  coefficients  can  be 
obtained  using  an  ionospheric  model  and  integrating  a  coupled  .set  of  differentitil  equations 
(equivalent  to  Maxwell's  equations)  for  the  field  components  E^,  Ey,  Hy  downward  from  a 
great  height  [Budden,  1961  j.  The  refiection  coefficient.s  are  funned  from  suitable  combinations 
of  these  field  components  below  the  ionosphere.  The  reflection  coefficients  depend  upon  the 
ionospheric  electron  density  and  collision  frequency  profiles,  the  geomagnetic  field  direction  and 
suength  as  well  as  the  ray  incidence  angle  and  wave  frequency.  Deteiinination  of  reflection 
coefficients  in  this  manner  is  called  a  "full  wave"  calculation. 


•  •  • 
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Figure  6-1 .  Schematic  of  the  LF  wave  hop  method  ray  path  showing  up-  and  down-going  rays  frohi  a 
lightning  source  at  -6  Km.  in  this  example  the  propagation  path  crosses  the  day-night 
terminator.  At  LF,  D-layer  reflections  occur  at  nominal  heights  of  70  and  90  km  for  day  and 
night,  respectively. 
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In  a  full  wave  calculation  four  complex  reflection  coefficients  (Rji,  R12,  R21,  R22)  obtained 
for  each  particular  set  of  iono.spheric  and  gei)magnetic  conditions,  wave  frequency,  and 
incidence  angle,  Figure  0-2  ,hows  typic.al  nighttime  reflection  coefficients  as  a  function  of 
frequency.  Rij  is  the  ratio  of  a  reflected  plane  wave  of  polarization  j  to  an  incident  plane  wave 
of  polarization  i  m.ea.sured  below  a  planar,  horizontally  stratified  magnetoionic  medium;  index  1 
denotes  verdcal  polarization  (magnetic  vector  transverse  to  plane  of  propagation),  index  2 
denotes  horizontal  polarization  (electric  vector  transver.se  to  plane  of  propagation).  The  four  Rij 


determine  a  matrix  R  = 


P,.  R 

R„  R 


12 


22 


.  For  multiple  hop  propagation,  the  overall  reflection  co¬ 


efficient  will  be  the  matrix  product  R*^  =  R^Cin  1  Rn  1  2>  •••’  RaGl^^l  where  the  matrices 

ro,,  0  1 

G  =  I  are  the  fresne!  reflection  coefficients  tor  the  intermediate  ground  reflection 

^0  ^22 

points.  In  particular,  the  magnitude  ||Rn||  of  the  1,1  component  of  the  product  mati'ix  RN  gives 


the  overall  path  loss  (due  to  reflections)  for  a  vertical  electric  dipole  transmitter  (VED)  and  a 
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Exact  refleciioii  coetiicients,  grazing  incidence 

Nighttime  ionosphere,  s,sn  =  0,  azimuth  to  mag.  field  =  45°,  dip  =  60° 
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Rgure  6-2.  Typical  nighttime  ionospheric  refiection  coefficients  computed  by  full  wave  integration 

[Sudden,  1961]  shown  as  a  function  of  wavo  frequency  for  grazing  ionospheric  incidence 
(about  80 )  and  fixed  magnetic  dip  and  declination  angles  (low  solar  activity). 

VED  receiver  .separated  by  N  ionospheric  reflections  and  N-1  intermediate  ground  reflections. 
Then  the  electric  field  of  thi;;  N  hop  mode  will  be 

=K.D-'G;,G;,.a^.F^.R,^  V/m  (fi.l) 

where 

K  is  a  proportionality  consuint  depending  on  the  U'ansmitter  power  in  watts  (tlie  ENT 
power) 

DfM  is  the  physical  length  of  the  N  hop  ray  connecting  source  and  receiver  in  km 

Gj,, ,  G are  the  transmitter  and  receiver  antenna  gain  patterns,  e.g.,  for  simple  vertical 
electric  dipole  =  cos  \\i  ^  where  1]/^  is  the  takeoff  angle  of  the  N  hop  ray. 

is  the  component  of  the  ov  all  N  hop  reflection  coefficient  matrix  for  source 
polarization  i  and  receiver  polarization  j. 


ttfsj  is  the  convergence  correcuon  discussed  in  Wait  [  1960]  and  JohI,r  [1961J  which 
converts  the  plane  wave  retlection  coet'ficieni  R*''  into  an  effective  spherical 
reflecdon  coefficient. 


ai^j  =  1^1  +  — j  2Nsin-^~/sin9  |  a|^l 


-cos-^l+h  /  {a  +  h)cos-^ — a 
2N;  j  L  2N 


A, H;;3'[Z^]exp  -i[— -z. 


h  =  reflection  height  (effective  ionospheric  reflection  height) 
a  =  radius  of  earth 

0  =  [transmitter-receiver  distance |/[ radius  of  earth) 


Zjsj  =  oyc  sin  VN/3cos2<j/|.j 
(i3  =  angular  wave  frequency 
c  =  speed  of  light 


\)/{M  =  N  liop  ray  takeoff  angle 

I-I|y3  =  is  the  iiankel  function  of  order  1/3,  of  the  second  kind  [  Abramovitz  and  Stegun, 
1972;  Cambridge  Computation!  Laboratory,  1945]. 


L'n  is  the  cutback  factor  discussed  in  Wait  and  Conda  1 1958]  which  accounts  for  the 
presence  of  the  earth  at  bo'h  transmitter  and  receiver.  In  the  geometrical  optics 
region  this  may  be  written 

I'N  =[i+R:(VN)][‘+f<;(VN)]  (6.2) 

Rfi  are  the  fresnel  plane  reflection  coefficients  of  the  ground  at  the  source  and  receiver, 
takeoff  angle  source  and  receiver  polarization  i  and  j.  A  ground  conductivity 
map  indicating  land  and  water  regions  of  the  earth  is  used. 

We  compute  the  electric  intensity  due  to  a  single  lightning  flash  as  the  RMS  sum 


M  ,  '« 

^RMS  ~  ^  I^nII  ^  h 


•  6 


•  « 


(over  some  sufficient  number  of  hop  models  m  -  presently  m  =  4)  for  both  up-  and  down-going 
rays  emanating  from  an  elevated  isotropic  point  source  (the  ENT). 

A  ground  wave  computation  is  also  performed  on  the  closest  lightning  sources  only  (local 
■storms)  and  the  squared  magnitude  of  the  ground  wave  is  summed  with  the  geometric  ray  modes 
above,  (see  Section  9.0  below) 

To  facilitate  practical  computations  we  constructed  a  data  base  of  precomputed  reflection 
coefficients  by  repeatedly  integrating  the  coupled  field  equations  [Budden,  1961]  through  model 
daytime  and  nighttime  ionospheres  (under  both  high  and  low  sunspot  conditions)  based  on  the 
IRI-90  ionosphere  in  Bilitza  et  al.,  [1990].  Figure  6-3  shows  model  day  and  nighttime  D  and  E 
layers  from  1R-I90.  The  computations  were  performed  for  10  wave  frequencies  between  60  and 
900  kHz;  10  incidence  angles  between  1  and  80°:  8  values  of  magnetic  field  azimuth  (with 
respect  to  propagation  direction)  between  -90  and  225°;  7  values  of  geomagnetic  field  co-dip 
angle  between  0  and  180°  (the  electron  gyrofrequency  was  fixed  at  1.3  MHz).  The  resulting 
data  base  serves  as  tlie  skeleton  for  a  four-dimensional  linr  -ir  interpolation  procedure  to  obtain 
reflection  coefficients  for  day,  night,  high  or  low  sunspot  activity  for  any  combination  of 
frequency,  incidence  angle,  magnetic  azimuth,  propagation  direction,  tmd  co-dip  angle.  Figure 
6-4  compares  the  exact  and  interpolated  reflection  coefficient  as  a  function  of  frequency. 


Height  (km) 


Figure  6-3,  Typical  day  and  nighttime  ionospheric  D  and  E-!ayer  profiles  from  the  IRI-90  model 

ionosphere.  Profiles  shown  are  for  a  variety  of  geographic  latitudes  and  are  averaged 
over  season. 


Figure  6-4.  Comparison  between  full  wave  daytime  ionospheric  reflection  coefficients  and  coefficients 
interpolated  from  pre-computed  database  (see  text).  Also  shown  is  the  empirical  data  from 
CCIR  Report  265-7  [CCIR,  1990], 

The  CCIR  Report  265-7  [CCIR,  1990|  contains  graphs  of  empirical  values  of  the  reflection 
coefficient  Rj  i  (vertical  source  and  receiver  polarization)  as  a  function  of  wave  frequency  and 
ionospheric  ray  incidence  angle  (the  functional  dependence  is  on  f .  cos(i)).  Curves  are  given 
for  all  seasons,  nighttime  ionosphere,  and  a  separate  summer  and  winter  daytime  set.  Also  given 
are  curves  which  indicate  the  variation  in  reflection  coefficient  magnitude  from  solar  minimum 
to  solar  maximum  conditions.  These  data  shown  in  Figure  6-5  were  based  on  numerous 
propagation  measuremcnrs  in  which  oblique  incidence  ionospheric  reflection  losses  were 
inferred  from  skywave-groundwave  interference.  The  CCIR  reflection  coefficients,  being 
averages  over  many  propagation  paths,  do  not  have  explicit  dependence  on  magnetic  field 
direction.  We  found,  however,  that  tlrcre  was  good  agreement  between  the  precomputed 
fullwave  rellection  coefficients  and  the  CCIR  curves  which  we  digitized  and  implemented  in  a 
computer  program  especially  for  the  case  of  grazing  incidence.  Figure  6-6  shows  a  typical 
comparison. 

From  the  standpoint  of  computational  efficiency  in  the  NoiseProp  program,  utilizing  the  CCIR 
reflection  coefficients  instead  of  the  fullwave  coefficients  (which  still  require  a  4D  interpolation 
and  magnetic  field  evaluation  at  every  ionospheric  reflection  area)  provides  a  significant 
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computation  time  savings.  Since  the  NoiseProp  program  sums  noise  propagated  over  hundreds 
of  different  paths  to  the  receiver,  it  appears  that  the  averaged-out  aspect  of  the  CCIR  reflection 
coefficients  has  a  minimal  effect  on  the  overall  combination  of  path  losses.  We  have  retained 
both  the  fullwave  and  CCIR.  reflection  coefficients  but  the  pre.sent  implementation  of  NoiseProp 
has  been  calibrated  using  the  faster  CCIR  reflection  coefficient  calculation  in  the  LF  propagation 
section. 

7.0  MF  PROPAGATION  METHOD  (-900  kHz  to  2  MHz) 

At  frequencies  above  ~900  kHz  and  low  takeoff  angles  (con'esponding  to  -2000  km  hop 
lengths)  the  CCIR  265-7  reflection  coefficients  no  longer  are  applicable.  Moreover,  the 
complicated  effects  inherent  in  MF  propagation  (the  sensitive  role  of  the  geomagnetic  field 
direction  and  magnitude;  and  significant  short-term,  day  to  day,  diurnal,  seasonal,  geographical, 
and  solar  cycle  variations  in  ionospheiic  loss)  make  a  simple,  physical  model  (such  as  the  wave 
hop  method)  underdetennined  in  both  the  number  and  scope  of  the  contributing  physical 
phenomena  which  can  be  modeled. 

Various  techniques  have  been  developed  to  model  MF  propagation  [PoKempner,  1980].  We 
have  cho,sen  a  relatively  simple  empirical  method  published  by  the  CCIR  as  Recommendation 
435-7  [CCIR,  1992b].  Method  435-7  was  based  on  statistical  analyses  of  field  strength 
measurements  for  266  propagation  paths  distributed  worldwide;  as  well  as  analysis  of  separate 
areas  for  which  propagation  paths  were  unavailable.  This  type  of  method  is  usually  derived 
from  a  multiple  regression  analysis  of  monthly  median  field  strength  measurements  from  which 
the  dependencies  on  variables  such  as  frequency,  path  length,  solar  activity,  geomagnetic  field 
vector  at  intrapath  control  points,  solar  zenith  angle,  etc.  are  determined.  The  predicted  skywave 
field  .'^'trength  (vertical  polarization)  is  given  by: 

Er=  V  f  G,  -L„  -t- A-201og|„p  -l()-’KKP-L,  (7.1) 


Jr: 


m 


•  • 


where 

E  is  the  annual  median  of  half-houriy  median  field  streng^'"’  in  dB  (q.V/m) 

V  =  P  f  G,  where  P  is  the  radiated  source  power  in  dB/(l  KW)  and  G  accounts  for 
transmitting  antenna  ele  vation  and  azimuth  gain  and  efficiency  (dB) 
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Gg  is  the  sea  gain  correction  important  only  when  the  transmitter  or  receiver  is  near 
seawater 

Lp  is  the  polarization  coupling  loss  (dB) 

Lp  sl80[36  +  e'  +  -2 

where  9  is  the  path  azimuth  in  degrees  from  the  geomagnetic  E-W  direction, 
j9|  <  90”,  I  is  tlie  geomagnetic  dip  angle,  Lp  is  evaluated  (using  a  13  coefficient 

expansion  for  magnetic  dip  and  declination)  at  both  source  and  receiver  locations  and 
the  results  added. 

A  106.6  -  2sin<I>  where  <&  is  a  path  averaged  value  of  geomagnetic/latitude. 

P  is  the  slant  path  length  for  E-layer  reflections  approximately  (d^  +  40,000) 
where  d  is  the  path  great  circle  distance  in  km. 

Kj^  is  the  loss  factor,  =  3.2  +  0. 19f  '^  tan^Cd)  -h  3)  +  0.01  •  b  •  SSN 

where  f  is  the  frequency  in  kHz,  SSN  is  the  12  month  running  average  sunspot 
number,  b  =  4  for  North  American  paths,  1  for  European  and  Australian  paths,  and 
0  elsewhere. 

Lf  is  the  diurnal  loss  factor. 

The  method  depends  on  graphs  contained  in  CCIR  report  435-7  [CCIR,  1992b]  for  the  sea  gain 
Gs  and  diurnal  loss  term  Lt.  Furthermore,  rather  involved  averaging  rules  are  given  for  the 
determination  of  a  single  value  of  the  geomagnei  ic  latitude  O  to  use  for  patlis  up  to  12000  km. 
The  algorithm  was  computerized  by  Inuki  et  ai.  1 1983]  and  validated  against  observations.  It  is, 
of  course,  also  intrinsically  consistent  with  the  CCIR  path  data  base  used  to  derive  it.  We  have 
programmed  a  slightly  simplified  version  of  the  algoritliin,  similar  the  summary  given  in  Davies 
[  1990,  pp.  433-436].  We  u,se  the  algorithm  with  its  weak  fy,j  loss  ilependencc  up  to  2000  kHz. 
Above  2000  kHz  an  HF  propagation  method  is  implemented.  As  in  the  LF  section,  the 
NoiseProp  algorithm  sums  MF  po  wer  radiated  by  using  250  lightning  ENT  sources. 


8.0  NOISEPROP  HF  PROPAGATION  METHOD 

HF  propagation  prediction.s,  unlike  LF  and  MF,  require  a  comparatively  detailed  ionospheric 
description.  Above  2  MHz,  the  NoiseProp  program  performs  ionospheric  propagation 


calculations  using  empirical  E  and  F2  layer  numerical  maps  in  combination  with  algorithms 
describing  the  takeoff  angle  and  virtual  reflection  height  of  ray  paths  from  source  to  receiver. 
The  propagation  algorithm  is  based  on  the  CCIR  Recommendation  533-3  [CCIR,  1992c]  method 
though  there  are  certain  significant  differences  in  our  approach.  'ITie  description  of  the 
ionospheric  F2  layer  parameters,  foF2,  M(3000)F2,  hmF2  are  those  used  in  the  International 
Reference  Ionosphere  IRi-90  [Bilitza,  1990]  and  the  CCIR  Recommendation  434-5  [CCIR, 
1992a].  This  F2  layer  description  is  tlie  so-called  CCIR  "Oslo"  ionosphere  [CCIR,  1966,  1967] 
represented  by  12  monthly  15  coefficient  fourier  series  diurnal  variation  Mid  a  spherical 
harmonic  expansion  in  terms  of  the  modified  dip  angle,  arctan  [l  /  yfeos  <j)] ,  where  I  is  the 

geomagnetic  dip  angle  and  ([)  is  the  geographic  latitude.  The  model  is  linearly  interpolated  for 
sunspot  numbers  between  0  and  100.  The  E  layer  critieal  frequency,  foE,  is  also  based  on  the 
formulation  used  in  IRI-90  and  CCIR  434-5  [CCIR,  1992a]. 

The  propagation  model  then  determines  several  ray  paths  connecting  the  source  to  the  receiver 
as  follows: 

•  The  great  circle  path  (short  path)  between  the  ENT  source  and  receiver  is  determined  and 
divided  into  equal  subsegments,  one  per  hop.  The  maximal  hop  length  is  determined  with 
respect  to  the  F2  layer  maximum  height,  hml"2  at  midpath  and  will  usually  be 
approximately  4000  km. 

•  For  each  sub.segment  a  one  hop  mode  is  sought  as  follows: 

—  An  F2  mode  is  first  sought  by  utilizing  the  foF2,  FoE,  hmF2  data  corresponding  to  the 
given  universal  time,  month,  and  sunsf  oi  number  at  the  segment  midpoint  (or  averaged 
over  die  midpoint  values  and  values  at  the  adjacent  points  a  quarter  hop  on  either  side). 

—  An  equivalent  triangular  reflection  height  (equivalent  virtual  height)  is  determined  using 
the  empirical  algorithm  described  in  CCIR  533-3.  'Ihe  empirical  algorithm  is  used  in 
place  of  the  slower  iterative  procedure  which  solves  the  transmission  equations  for  a 
double  layered  parabolic  ionosphere  to  detemrine  ray  takeoff  angle  and  virtual  height  as 
a  function  of  path  length  [CCIR,  19H6|.  The  algorithm  detennines  virtual  height  as  a 
fund  fin  of  frequency,  hop  length,  and  the  state  of  the  ionosphere  (E  and  F  layers)  along 
I  ho  hup  segment. 

•  The  elevation  angle  thus  determined  for  an  F2  hop  is  te.sted  for  viability.  In  this  ca.se,  we 
check  to  see  if  the  takeoff  angle  is  positive.  If  not,  the  hop  segment  is  subdivided  and  the 


process  repeated  until  an  F2  mode  of  positive  takeoff  angle  is  found  for  each  path 
subsegment. 

•  The  F2  mode  on  each  subsegment  is  tested  for  E  layer  reflection.  If  tlie  E  layer  critical 
frequency  at  the  midpoint  of  any  path  subsegment  is  sufficiently  high,  foE  >  f/1.05  sectj), 
where  f  is  the  wave  frequency  and  (j)  is  the  ray  incidence  angle  at  1 10  km.  The  F2  mode  is 
completely  occulted  by  tlie  E  layer,  rn  this  case  the  F2  hop  on  this  subsegment  is  replaced 
by  one  or  more  E  hops  with  virtual  reflection  height  at  1 10  km. 

•  At  this  stage  the  path  from  the  source  to  receiver  has  been  subdivided  into  E  and  F  layer 
hops.  The  E  and  F  layer  MUF  (maximum  usable  frequency)  are  determined  using  the 
empirical  methods  in  CCIR  434-5  [CCIR,  1992al  and  CCIR  533-3  [CCIR,  1992cJ.  These 
are  used  to  formulate  a  circuit  MUF.  In  communications  prediction,  the  circuit  MUF  is  used 
to  determine  a  circuit  availability  percentage.  In  our  NoiseProp  application,  the  circuit  MUF 
is  used  to  determine  to  determine  an  "above  the  MUF  loss"  term  used  in  the  subsequent 
propagation  loss  budget. 

•  The  pass  ranges  and  virtual  reflection  heights  are  sent  to  an  adjustment  algorithm.  Here,  the 

virtual  heights  are  fixed,  while  the  pass  ranges  are  adjusted  so  that  incident  and  reflected  ray 
angles  are  equal  at  each  ground  reflection  point.  This  adjustment  forces  the  complete  mode 
to  be  geometrically  consistent  and  provides  an  adequate  approximation  to  an  exact  physical 
solution.  The  exact  solution  relating  the  takeoff  angle  of  a  ray  connecting  a  receiver  and 
transmitter  (separated  by  multiple  hops  .spanning  varying  ionospheric  conditions)  to  the  total 
distance  spanned  can  only  be  found  by  a  time  consuming  iterative  ray  tracing  search 
(shooting).  In  principle,  the  udjuste  '  pass  ranges  could  be  returned  to  the  virtual  height 
algorithm,  new  virtual  heights  and  'akcoff  ant'les  determined,  and  the  resulting  ranges  and 
ranges  sent  to  the  adjustmen.  al  jo;  '  in  a  uniform  takeoff  angle  close  to  the  true 

solution.  This  iteration  would  i  oir  ,.gc  U)  the  exact  physical  solution.  The  accuracy  of  the 
first  adjustment  is  sufficient,  in  ''.  aver,  for  all  communications  oriented  predictions  and 
therefore  the  Noiseftop  ENT  r  upagation  model  as  well.  Tire  adju.stnient  is  illustrated  in 
Figure  8- 1 . 

The  median  RMS  field  strength  from  a  single  ENT  source  in  dB  above  1  |xV/m  is  the  root-sum- 

squaie  of  individual  modes; 

Etot  (8.1) 

k=l 
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a.  Before  adjustment 


b.  After  adjustment 


Figure  Q-1 .  Pass-range  adjustment  in  HF  virtual  raylrace  propagation  algorithm.  Virtual  heights  are 
unchanged  from  a.)  to  b.). 
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where  is  the  median  RMS  field  .strength  of  the  K'**  muitihop  mode  (in  dB/ljiV/m)  and 
Ej.  =  1 07.2 -f  20 log, „  f(MHz)  +  P,,;^,,  H-G(a)-L(a)  (S.2) 


where 

f  =  wave  frequency  in  MHz 
Pent  -  source  power  in  dB  above  I  watt 
a  =  the  adjusted  elevation  angle  of  the  mode 
G(a)  =  source  gain  at  elevation  angle  a  (dB) 

L(a)  =  propagation  los.s  of  K^’  mode 

I’re.sently,  the  ENT  source  is  modeled  as  an  isotropic  radiator  and  so  G(a)  =  0.  L(a)  accounts 
for  all  losses  accrued  by  the  mode  along  its  path.  This  term  is  composed  of  the  following  parts; 


'^  ^abs  ■^Snd  '^^pol  ^hf  ^  ^Ip '^^MUF 


(8.3) 


•  L(,j-  is  the  basic  free-space  loss  due  to  the  inverse  distance  falloff  of  the  electric  field 

magnitude.  L|,f  =  20  logioKTtp'A,]  where  p'  is  the  ray  path  length  and  X  is  the  wavelength. 
In  kilometer  and  mflz  units:  =  32.45  +  20  •  log  f  +  20  log  p'. 

•  is  the  non-deviative  absorption  loss  due  to  D  and  E  layer  tiansits.  There  are  various 
empirical  models  which  reduce  the  total  jibsorption  so  a  sum  of  terms,  one  per  hop,  or  one 
per  control  point  e.g.,  CCIR  533-3  [CCIR,  1992c],  Barghausen  et  al.  f  1969],  Lloyd  et  al. 
[1978],  Pricker  [1987],  and  Sailors  and  Rose  [1993].  All  die  empirical  method.s  except  for 
that  developed  by  the  BBC  [Pricker,  1987]  are  control  point  methods  and  average  absorption 
only  over  4  to  5  control  point  locations.  Our  virtual  ray  trace,  however,  depends  explicitly 
on  the  ionospheric  conditions  from  hop  to  hop.  We  therefore  u.se  the  BBC  absorption 
formula  which  e.\plicitly  sums  the  ionospheric  absorption  from  each  leg  of  each  hop 


L 


eabi 


5]^13-»<D 


(foE)^  scci|, 
(f  +  ().6)^ 


(8.4) 


where  the  sum  is  over  up  and  downgoing  legs  of  all  hops,  foE  is  the  E-Iayer  critical 
frequency  at  the  point  of  inter.section  of  the  ray  and  1 15  km  altitude,  ig  is  the  incidence  angle 
of  the  unrcfracted  ray  at  1 15  km  (|i,,  s  cos"'  (l  -  0.965 cos'  j  and  takes  into  account 

the  depth  of  penetration  of  the  ray.  For  an  P2  hop 


d),.  =  1  +  • 


0.07 


foE  sec  i  I 


•0.75 


(8.5) 


For  an  E-layer  hop 


1.2 


0.41  . 


(8.6) 


It  is  assumed  that  E-layer  reflection  occurs  when 

- - <0.95  . 

foE  .sec  i ,, 
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•  Lgjuj  is  the  loss  due  to  ground  reflection  on  multi-hop  modes.  A  1“  x  1°  land-water  map  is 
used  to  determine  the  relative  dielectric  con.stant  and  conductivity  o  at  each  ground 
reflection  point  (seawater:  =  80,  o  =  5;  land:  =  0.4,  o  =  0.00 1 ).  The  Fresnel  reflection 
coefficients  for  vertical  and  horizontal  polarization  are  computed 

R„  =  sina- -  cos^  a)'^^j/j^n’ sina  +  (n* -cos^  (8.7) 

R|,  =^sina-(n^  -cos^  a)'^^j/j^sina4-(n^ -cos^  a)'^^  (8.8) 


where  =  Ej.  -  i(18000CT)/f  MHz  is  the  complex  refractive  index  of  the  ground,  a  =  ray 
takeoff  angle.  Then,  assuming  random  polanzation  of  the  downgoing  wave  (a  valid 
assumption  for  a  ray  after  one  or  two  ionospheric  transits)  we  set  the  ground  loss 


k..a=10.1ogJ^ 


Kf + Rjr 


Lpol  is  the  polarization  coupling  los.s  due  to  the  iiiiperfect  coupling  between  the  initial  wave 
polarization  and  the  two  chaiacteristic  polaiizations  preferred  by  the  ionosphere  (X  and  O). 
Wo  use  a  simple  empirical  formula  described  in  Sailors  and  Ro.se  1 1993|: 


Lpoi  =  A  4-  B  .  (average  hop  length  in  radians)  -f  Nhops 


where 


A  =  (2.400.  Ne  +  2.060  .  N^,-)/N  hops 
B  =  (2.548  .  Nh  +  1.019  .  Np)/N  hops 
Np  =  number  of  E  hops 
Np  =  number  of  F  hops. 

N  hops  =  Np  -r  Np 

is  the  excess  system  loss  from  the  tables  in  Barghausen  et  al.  iT9u9j,  CCIR  252-2  and 
533-3  [CCIR,  1986;  1992cl.  This  empirical  term  allows  for  auroral  and  other  high  latitude 
losses.  It  is  evaluated  in  terms  of  iniilpath  geometric  latitude  (N  or  S  of  equator),  season,  and 
midpath  local  time  and  is  0  for  geometric  latitudes  ie,ss  than  42.5°  (also  see  Lloyd  et  al., 
1978). 
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Lmuf  is  the  "above  the  MUF"  loss.  The  model  does  not  assume  that  propagation  ceases 
above  the  geometrical  optics  path  MUF  (junction  frequency).  It  is  known  that  the  signal 

docs  not  fall  completely  to  zero  at  frequencies  above  the  MUF.  CCIR  533-3  [CCIR,  1992el 
gives  the  following 


0  if  f  <  F  MUF 


for  E  modes 


Lmui.-  -  min  81dB,13()* 


EMUF 


(8.10) 


for  F  modes 


0  if  f  <FMUF 


min  62dD,36«  — - - --1 

FMUF 


(8.11) 


The.se  emi)irical  curves  account  for  both  1 )  tlie  effect  on  the  median  monthly  field  strength 
due  to  the  source  frequency  exceeding  the  daily  MUF  on  some  day.s  and  not  propagating,  and 
2)  losses  occurring  (due  to  phenomena  other  than  simple  two-dimensional  ray  optics)  on 
days  when  the  source  frequency  i.s  above  the  basic  MUF  but  propagation  occurs  nonetheless 
(due  to  non-classical  or  non-ray  optics  phenomena). 

Gijj'is  a  horizon  focus  gain  for  short  paths  (less  than  1/4  of  the  earth’s  circumference, 

- 1(1,(100  km).  It  is  a  polynomia!  fit  (in  elevation  angle)  to  curves  ba.sed  on  double  parabolic 
E,  F  layer  ionospheres.  It  leaches  a  maximum  of  9  dU  at  small  takeoff  angles  as  shov/ii  in 
Figure  8-2.  The  curves  are  given  by  Sailors  and  Rose  1 1993|: 

6 

G|,,  =  J^ajX'  for  elevation  angles  <  l(F,  x  =  (elevation  angle  -5")/5'^, 

ii£0 

ft 

G,,,  =  ]^biy'  for  elevation  angles  >10",  y  =  (elevation  angle  -50")/40”, 

i=0 

and  the  coefficients  are  given  by 
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lW»»tT*Tn'4«Til 


'‘H) 

ai 

a2 

“3 

^^4 

as 

r 

E- modes 

4.14 

-4.00 

3.007 

0.069 

-1.06 

0.681 

-0.347 

F-niodcs 

6.03 

-3.161 

1.400 

0.624 

-1.413 

0.088 

0.533 

F-modes  I  1.50 


•0.274  0.387 


0.540  -0.267  -1.136  0.853 


Gip  is  a  long  path  I'oous  gain  which  is  effective  for  circuit  lengths  greater  than  1/4  of  the 
earth's  circumference  and  i.s  given  by 


G|,.  =20  logio  l- 


6371k 

D 
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where  D  is  the  circuit  path  length  in  km  (always  less  than  one-half  the  earth's  circumference) 
[Hortenbach  and  Rogler,  1979J. 

8.1  Influence  of  Geomagnetic  Activity 

Extensive  studies  have  shown  that  the  F2  layer  maximum  critical  frequency  foF2  and  the  F2 
layer  MUF  are  significantly  correlated  with  measured  values  of  the  local  geomagnetic  k-index 
[Davis,  1970;  Wrenn  and  Rodger,  1989].  Two  tendencies  have  been  noted  during  periods  of 
increased  geomagnetic  activity  1)  a  relative  reduction  of  foF2  at  auroral  latitudes  in  summer  or 
during  high  sunspot  activity  in  local  daytime,  and  2)  a  relative  enhancement  of  foF2  at  low 
geomagnetic  latitudes  (above  35°)  in  winter  or  during  low  sunspot  activity  in  the  local  afternoon 
and  night.  At  iniddle  latitudes  these  two  tendencies  are  combined  to  a  degree  depending  on  tiie 
particular  geomagnetic  latitude,  local  time,  and  season. 

Barghausen  et  al.,  [  1969]  used  the  results  of  Davis  [  1970]  to  adjust  a  given  path  MUF  by  a 
fraction  P  =  (actual  MUF)/(montlily  median  MUF)  for  a  given  local  time  and  geomagnetic 
activity  index.  7’he  fraction  F  is  given  by  a  linear  regression  on  geomagnetic  activity  index  k: 

P  =  P„  +  bk  (8.12) 

where  Po  and  b  are  coefficients  which  depend  on  geomagnetic  latitude  at  midpath,  tluce  levels 
of  sunspot  number,  three  seasons,  and  eight  local  time  blocks.  The  coefficients  are  linearly 
interpolated  in  sunspot  number.  The  coefficient  of  correlation  c  from  the  regression  analysis  is 
also  stored  and  the  adjustment  Actual  MUF'  =  P  •  (monthly  median  MUF)  is  only  made  for 
conditions  where  lei  >  0.3. 

The  adjustment  of  the  path  MUF  for  k  index  will  affect  the  abovc-the-MUF  loss  term  in 
the  propagation  loss  budget,  particularly  for  propagation  paths  which  cross  high  geomagnetic 
latitudes.  In  tlie  NoiseProp  summation  over  multiple  ENT  .sources,  the  effect  of  the  adjustment 
is  probably  only  seen  at  high  latitudes  when  the  frequency  is  margin-ilIy  greater  than  the  path 
MUF  for  a  majority  of  paths  between  the  dominant  ENT  soutces  and  the  receiver.  The  k  index 
adjustment  can  be  defeated  altogether  in  the  NoiseProp  inputs:  in  which  case  the  MUFs  and  the 
corresponding  loss  terms  will  all  reprc.sent  median  conditions. 


9.0  GROUNDWAVE  PROPAGATION  OF  LOCAL  ENTs 


For  local  lightning  ENT  sources  (closer  than  ~5()()  km  to  the  receiver)  we  must  consider 
propagation  by  the  groundwave.  The  groundwave  field  derives  from  the  solution  to  tlie  wave 
equation  over  a  smooth,  finitely  conducting  spherical  earth.  The  source  will  be  a  vertical  electric 
dipole  at  an  elevated  point  height  h  above  the  surface.  The  VED  is  the  most  efficient  exciter  of 
the  vertically  polarized  groundwave  .  (An  HED  can  also  excite  a  groundwave  but  it  will  be 
orders  of  magnitude  weaker  and  irrelevant  from  the  standpoint  of  our  NoiseProp  Mimmation  of 
ENT  power.)  Groundwave  propagation  is  most  efficient  at  lower  frequencies,  shorter  distances, 
and  higher  ground  conductivities. 

The  vertical  electric  field  at  great  circle  distance  d  =  a  •  0  measured  along  the  earth's  surfice 
[Wait,  1970,  1974;  Bremmer,  1949J  is 

E  =  Eo.W  (9.1) 

where  £„  is  the  unattenuated  field  of  a  VED  over  a  flat  perfectly  conducting  plane  and  W  is  the 
groundwave  attenuation  function.  ITiese  may  be  written 

lE„|=--^VPjv/m  (9.2) 
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E|.  =  relative  dielectric  C(5nstant  (ground) 
a  =  ground  conductivity 
i  =: 

tg  =  ^^(q)  are  roots  of  the  differential  equation 

dt.  1 


witli  the  initial  condition 
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where  the  otj;  satisfy  Aifa^)  =  0  s=l,2,3i.-,  and  Ai(z)  is  the  standard  Airy  function. 


Wi(t.  -  y,) 


and 


i(t.  -y^) 

w.(0 


are  the  height  gains  for  the  source  and  receiver.  Since  the 


receiver  is  on  the  ground,  y2  =  0  and  wc  replace  the  second  height  gain  factor  by  1  and 
Wi(t)  =  V7c[Bi(t)  -  i(Ai(t)]  where  Ai(t)  and  Bi(t)  are  the  standard  airy  functions  1  Abramovitz 

and  Stegun,  I972j. 

Tlte  expression  for  the  spherical  earth  attenuation  function  W  reduces  in  the  case  hj  =  h2  =  0  and 
ti/a  «  1  to  the  Somnierfeld  flat  earth  attenuation  function  F(p)  where 

b'(p)  =  1  -  i(7cp)*^^  c‘1’  erfefip*/^)  and  p  =  -i(kd/2)A2  (9.4) 

Groundwave  propagation  carves  for  a  1  kw  source  at  grouiul  level  for  frequencies  from  15  kl  Iz 
to  5  MHz  have  been  plotted  in  b'igure  9-1  for  coinpaiisoii  with  the  unatteiiuated  perlectly 
conducting  ground  dipole  pri)pagation  (top  curve).  I’he  higher  frequencies  are  strongly 
attenuated  beyond  a  few  hundred  kilometers. 

'Hie  NoiseProp  system  only  utilizes  die  groundwave  algorithm  for  local  sources  wliich  we  define 

kvj  uc>  wiuiiii  .'juiiiu  iK)  /  l\/  k^V9i.ii.Oui  Oi  iturkiiiaic  ij>auii  duuCuii  iti 

a  flashrate  density  in  accordance  with  that  weighting  scheme.  The  BNT  radiated  power  derived 
from  these  5°  x  5°  sources  is  Uien  propagated  by  the  groundwave  formalism  to  the  receiver 
(using  the  geometric  mean  ground  conductivity  and  dielectric  constant  along  die  patli)  and  tlie 
received  powers  summed. 
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propagated  by  the  groundwave  formalism  to  the  receiver  (using  the  geometric  mean  ground 
conductivity  and  dielectric  constant  along  the  path)  and  the  received  powers  summed. 
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FIguro  9-1 .  Ground  wave  signal  strength  In  micro  volts-per-ineter  over  average  (spherical)  earth  for 
a  1  kilowatt  vertical  electric  dipole  (VED)  at  ground  level.  Dashed  curve  shows  1/D 
falloff  for  1  kW  dIpolo  over  a  perfectly  conducting  plane  surface. 

10.0  NOISEPROP  PROGRAM  AND  POWER  SUMMATION 

10.1  Summary 

Given  a  frequency  (between  fit)  kHz  and  30  MHz),  a  date  (year,  montli,  day),  universal  time, 
sunspot  number,  and  receiver  location,  the  NoiseProp  program  determines  the  noise  factor  Fa  in 
decibels  relative  to  kTb  (k-'Boltzmann's  constant,  T  is  288‘’K,  b  is  eflecti'  e  receiver  noise 
bandwidth  in  Hz)  in  a  short  vertical  montrpole  over  a  conducting  ground  plane.  This  .igure  can 
be  computed  for  a  percentage  of  time  between  ().{)()  1  and  99.999%. 

The  NoiseProp  program  operates  in  a  scries  of  steps  illustrated  in  Figure  10-1. 
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1.  Tlie  appropriate  ■  easonal  and  diurnal  GLO  maps  (or  DGLO  maps  if  available)  are  converted 
into  corresponding  source  power  maps  or  ENT  maps. 

2.  A  sorting  algorithm  arranges  the  globally  distributed  ENT  maps  in  a  rough  dominance 
ordering. 

3.  Depending  on  frequency,  the  appropriate  propagation  algorithm  (LF,  ME,  HE)  is  applied 
repeatedly  to  each  ENT  source  in  order  of  dominance.  'I’he  available  power  in  the  vertical 
electric  field  at  the  receiver  location  is  summed  over  a  large  (presently  250)  number  of 
sources. 

4.  lire  local  ENT  power  in  the  nearest  four  5°  x  5°  subcells  is  propagated  by  the  groundwave 
algorithm  and  summed  at  the  receiver  location. 

5.  The  summed  available  power  at  the  receiver  is  converted  to  tire  decibel  value  Fa.  Fa  is 
treated  as  the  median  of  the  nortnal  distribution  from  which  the  level  exceeded  a  chosen 
percentage  of  time  is  computed. 
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Figure  10-1 .  Flow  chart  for  NoisoProp  computations. 


10.2  Discussion 


1.  ITie  conversion  of  GLO  or  DGLO  map  flashrate  data  into  ENT  map  power  data  is  a 
weighted  sum  of  return  stroke  and  bipolar  in-cloud  stroke  spectra  for  frequencies  below 
1  MHz.  Above  1  MHz  a  complete  flash  spectrum  is  utilized  as  discussed  in  Section  4.0. 

Flash  power  standard  deviations  are  calculated  as  in  Warber  and  Field  [1991].  This  requires 
weighted  sums  of  stroke  power  variances.  Above  1  MHz  we  do  not  use  separate  stroke 
phenomenology  for  the  u  .'erall  flash  power.  Due  tc  the  lack  of  data  on  flash  power  standard 
deviations  (at  HF)  comparable  to  that  available  at  VLF  (and  used  in  LNP)  we  have  carried 
over  the  LNP  formalism  to  compute  standard  deviations  above  1  MHz.  Tlie  overall  flash 
spectrum  (straight  line  in  Figure  4-1)  is  converted  by  fractional  multipliers  into  a  sum  of 
fictitious  .stroke  spectra  (return  stroke,  horizontal  and  vertical  bipolar  in-cloud)  for  which 
standard  deviation  expressions  are  available.  Then  the  standard  deviation  can  be  built  up  for 
the  overall  flash  using  the  existing  LNP  formalism.  The  overall  flash  variances  from  each 
ENT  are  summed  together  to  give  a  total  noise  variance.  The  total  variance  is  used  to 
determine  the  distribution  (assumed  log-normal)  of  the  overall  noise  power  and  the 
corresponding  exceedence  level  for  a  given  percentage  of  time. 

2.  The  preliminary  sorting  of  ENTs  into  a  dominance  order  is  done  by  comparirig  PgisiT/d^ 
where  Pent  ^NT  power  in  watts  and  d  is  the  distance  from  the  ENT  center  of  mass  to 
the  receiver  in  km.  This  very  ro uglily  orders  iNTs  according  to  thefr  propagated 
streng'ii  if  they  propagated  in  free  space.  We  found  it  recessary  to  add  up  a  large  number  of 
sources  from  this  presorted  orderL'.g  using  the  actual  LF,  MF,  or  HF  propagation  algorithms. 
This  is  because  the  propagation  functions  are,  in  fact,  not  monotonic  with  distance, 
especially  at  HF  where  the  .skip  distance  phenomenon  is  pronounced.  Originally,  we  had 
intended  to  simply  observe  the  size  of  the  power  sum  as  each  propagated  ENT  was  added 
and  to  temiinate  die  calculation  when  the  increment  in  the  sum  was  smaller  than  some  small 
fraction  of  1  dB.  However,  we  found  many  cases  where  although  the  stopping  criterion  was 
set,  subsequent  ENl's  in  the  dominance  ordering  contributed  significant  power  to  the 
cumulative  sum.  Therefore,  at  present,  we  add  the  first  250  (out  of  648)  ENTs  in  every 
calculation.  Although  this  number  i  >  sufficiently  lai  g.e  not  to  miss  any  significant  ENTs,  it  is 
admittedly  a  serious  limitation  in  the  time  economy  of  the  NoiseProp  program. 

3.  The  sky  wave  and  groundwave  algorithms  in  MoiseProp  all  compute  the  RMS  field  suength 
E  in  dB  (pV/m).  Then  e  =  10^/2^  is  the  RMS  field  strength  in  the  vicinity  of  the 
receiving  antenna  in  units  of  pV/m  and  the  associated  power  flux  is  p  =  in  watts/m^ 


where  Z(j  -  12071  is  the  impedance  of  free  space.  The  available  power  at  the  terminals  of  the 
receiving  antenna  is  then 

Pavail  =  P  • *  S  (lO'l) 


where  p  is  the  power  flux,  A  is  the  effective  isotropic  antenna  aperture  area  and  g  is  the 
receiving  antenna  gain  [Stutzman  and  Thiele,  1981].  Since  A  =  ^2/471;  where  X  is  the 
wavelength  in  meters,  we  have 


■“  Z„  *  47: 


4807t' 
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-(watts).  (10.2) 
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The  gain 's  nominally  set  to  unity  and  the  numbers  Pavail  summed  over  all  contributing 
ENT  sources  to  yield  the  overall  noise  power  in  watts  at  the  receiving  antenna  terminal. 
(Note:  if  we  express  dB  values  by  using  capital  letters,  the  above  equation  becomes 

Pavail=E  +  G-  107.2  -  20  logiof  MHz  O0.3) 

which  is  the  commonly  used  relation  seen  in  ail  the  standard  propagation  prediction  codes.) 


11.0  CALIBRATION  OF  NOISEPROP  MODELS 

The  NoiseProp  system  consists  of  two  sections;  the  ENT  source  algorithm  and  the  propagation 
algorithms  (LF,  M,  HE,  groundwavc).  The  propagation  algorithms  can  be  verified  to  a  certain 
extent  by  comparison  with  certain  existing  well  tested  propagation  codes,  hi  the  case  of  the  HF 
section,  extensive  comparisons  were  made  for  many  t:irc'’ii'>  y.'''tli  the  lONCAP  [Lloyd  et  al., 
1978]  and  ASAPS  [see  Sailois  lUid  Fo.se,  i993j  prediction  codes.  Tlie  codes  are  considered 
benchmarks  because  they  have  been  checked  and  adjusted  against  real  measurements  for  many 
years.  The  LF  wave  hop  code  can  be  tested  to  a  certain  extent  against  the  empirical  CCIR  435-7 
method  which  works  down  to  150  kHz. 

Tne  largest  source  of  uncertainty,  however,  is  the  lack  of  precise  data  on  lightning  flash  power 
spectra  above  the  VLF  range,  'lliis  introduces  a  systematic  error  into  all  results  which  is  a 
function  of  frequency  and  which  can  be  corrected  for,  a-posteriori. 
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A-posteriori  corrections  were  introduced  by  a  set  of  comparisons  with  the  CCIR  322-3  [CCIR, 
1988]  noise  model  at  the  50th  percentile  level.  The  CCIR  model  was  based  on  noise  data  taken 
at  16  stations  during  the  period  1957  through  1966  (not  all  the  stations  produced  data  for  the 
entire  period).  Tlie  data  were  grouped  into  four  .seasonal  blocks  and  six  four  hour  time  blocks 
during  each  day  of  the  seasonal  block.  The  CCIR  model  contains  median  hourly  values  of  RMS 
noise  power  in  each  time  block;  the  median  is  taken  v/ith  respect  to  the  seasonal  (~90  day) 
aggregates.  The  16  station  median  values  are  smoothly  inteipolated  onto  a  worldwide  grid  using 
spherical  polynomial  expansions  [Spaulding  and  Washburn,  1985].  The  coefficients  in  this 
expansion  are  the  basis  for  a  computerized  version  of  the  maps  in  CCIR  322-3. 

Systematic  conections  to  tlie  NoiseProp  predictions  take  the  form  of  a  function  of  both 
frequency  and  local  time  (at  the  receiver  site).  The  frequency  dependence  offsets  errors  in  the 
ENT  source  power  spectrum  (which  will  factor  systematically  through  all  subsequent  stages). 
The  local  time  correction  is  to  offset  any  cumulative  systematic  en'ors  in  the  propagation 
algorithms  tliem.selves  inasmuch  as  these  algorithms  all  contain  strong  local  time  variations  in 
their  propagation  loss  modeling.  A  local  time  correction  also  allows  for  accounting  of  non- 
specific  diurnal  trends  in  local  lightning  activity  not  explicitly  included  in  the  stroke  and  flash 
rate  models. 

The  correction  terms  were  computed  by  averaging  the  dB  difference  between  CCIR  322-3  and 
NoiseProp  over  four  montlr  ensembles  of  12  CCIR  station  measurements  in  4-hour  blocks.  Tire 
station  sites  ae  shown  in  Figure  11-1.  'Hie  output  of  the  CCIR  322-3  model  is  assumed  to  be 
closest  to  the  underlying  data  at  the  data  collection  location.s  themselves  and  at  greater  variance 
with  reality  in  areas  far  from  the  16  station  measurement  network,  e.g.,  over  oceans  [Sailors, 
1993],  Tlie  four  months  were  January,  April,  July,  and  October;  hence,  seasonal  effects  are 
averaged  out.  Since  the  averages  were  over  12  widely  .separated  locations;  station  .specific 
effects  are  also  averaged  out.  These  ensemble  averages  were  computed  at  100,  500,  1000,  and 
1500  kHz  to  calibrate  the  LF  and  MF  .sections;  5,  1(',  15.  and  20  MHz  to  calibrate  the  HF 
section.  The  resulting  calibration  curves  are  smoothly  approximated  by  least-.squares 
polynomials  in  local  time  and  linearly  interpolated  as  functions  of  frequency.  That  is,  we 
compute 

Adocaltime,  fkl-lz)-<NoiseProp-CClR>  (dB)  .  (11.1) 
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Tliese  calibration  curves  are  shown  for  the  eight  base  frequencies  in  Figures  1 1-2  and  1 1-3. 

The  correction  term  A  is  then  subtracted  from  the  uncalibrated  NoiseProp  output  after  each 
propagation  calculation.  Examples  of  NoiseProp  output  after  calibration  are  shown  in  Figures 
11-4  and  11-5.  Figure  11-4  is  a  set  of  LM-HF  noise  calculations  for  July,  0000  UT  and  low 


100  kHz  500  kHz  lOOCkHi  1600  kHz 


Figure  1 1  -2.  LF-MF  NoiseProp  calibration  curves  at  100,  500, 1 000,  and  1 500  kHz.  A  smooth 
polynomial  fit  is  superimposed  upon  the  actual  data.  These  curves  are  subtracted 
from  the  NoiseProp  output  to  account  for  systematic  frequency  dependent  and 
local-time  dependent  deviations.  The  curves  are  linearly  Interpolated  in  frequency. 
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Rgure  11-3.  HF  NoiseProp  calibration  curves  at  5, 1 0,  1 5,  and  20  MHz. 
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sunspot  activity.  The  solid  curves  are  the  Noi.seProp  calibrated  values;  the  dots  shown  CCIR 
values,  Figure  1 1-5  is  a  typical  set  for  HF  at  about  OXOO  UT.  The  graphs  show  the  result  of  100 
calculations  at  10  longitudes  from  0  to  3fi0°;  at  each  longitude  calculations  were  done  at  10 
latitudes  between  -60  and  60°;  this  sequence  of  c.tlculations  is  plotted  in  natural  ordering  of 


given  by  K  =  1,2,3,...,  100  where  K  =  10  •(i-1)  +J. 
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The  pattern  of  minima  at  tire  extremes  of  latitude  (±60°)  and  maxima  at  the  equator  is  seen  in 
botli  die  NoiseProp  and  CCIR  models.  Ihe  minima  in  the  CCIR  plots  are  not  as  deep  as  the 
NoiseProp  minima.  This  may  relicct  the  fact  that  the  CCIR  model  is  consti'ained  so  as  to 
smoothly  fit  a  sparse  set  of  data  (16  points)  to  the  worldwide  grid.  Ftir  from  data  points,  the 
properties  of  the  fit  ore  largely  determined  by  tire  nature  of  the  fitting  expansion,  which  is 
smooth  and  slowly  varying  radier  (liari  by  the  underlying  physical  phenomena  being 
approximated  ISailors,  10031. 

Ideally,  testing  of  a  source  based  noise  model  such  as  LNP  rrr  Noi.seProp  should  be  done  with 
source  distributions  representing  the  actual  distribution  of  lighlniug  sources  during  tire  period  rrf 
a  noise  measurement.  Cletu'ly  this  cannot  be  apirlied  to  the  CCIR  322-3  motlel  which  represents 
long  term  median  noise  levels  measured  over  a  great  variety  of  (many  years)  meteorological 
conditions.  Nor  is  it  at  all  completely  obvious  that  tire  median  noise  ineasurenrents  arc 
necessarily  congruent  to  a  NuiseProp  prediction  utilizing  median  (Itmg  tenir)  GLO  maps,  lire 
nonlinearity  in  the  propagation  meclrauism  will  alter  the  probability  di-stribution  of  ENT 
intensities  when  it  is  transformed  to  a  pi  obability  di.strihiition  of  noise  power  summed  over  many 
sourees,  'Hie  best  program  for  validation  and  calibration  requires  an  active  field  measurement 
program  simulUmeous  with  the  processing  of  current  meteorological  data  into  dynamic  IXJLO 
maps  which  are  input  to  Noi.seProp  (and  l.NP). 

12.0  NOISEPROP  APPLICATIONS  USING  DYNAMIC  LIGHTNING  OCCURRENCE 
MAPS 

'lliis  section  describes  the  development  of  an  operational  dynaniic  liglitning  forcc.asting  system. 
Such  a  system  will  produce  tlie  DGLO  maps  from  meteorological  data  which  can  serve  as  input 
data  to  the  LNP  or  NoiseProp  noise  prediction  codes.  'Ihc  replacement  (rf  the  long  term  median 
GLO  maps  (derived  from  satellite  obseiwations)  with  dynamic  DGLO  maps  will  be  tran.sparcnt 
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for  the  Operation  of  the  NoiseProp  code;  the  GLO  and  DGLO  map  files  will  have  identical  data 
fonnats.  The  DGLO  maps  will  contain  DGLO  forecast  based  data  in  those  sectors  where 
adequate  meteorological  data  is  available;  they  will  ettntain  the  original  GLO  data  in  sectors 
where  forecasts  could  not  be  made  (due  to  insufficient  data  at  a  given  time  and  place). 

We  have  developed  several  techniques  for  converting  short,  medium,  and  long  term  weathc,r 
forecasts  (along  with  other  parameters)  to  produce  lightning  flash  rates.  These  predictions  will 
be  assembled  into  the  DGLO  maps  and  be  made  tivailable  to  users  by  downloading  from  a 
central  location.  Only  the  military  prtKluccs  weather  forecasts  on  a  worldwide  basis.  It  appears 
tliat  access  to  data  from  the  Air  Force  Global  Weatlier  Cential  (AFGWC),  in  particular,  will  be 
needed  on  a  regular  basis  for  Itigh  confidence  forecasts.  A  detailed  feasibility  study  for  the 
lightning  forecasting  problem  is  given  in  a  companion  report  by  Wai'ber  and  Sinclair  1 1995]. 


12.1  Correlation  Between  Flash  Rate  and  Weather  Parameters 

In  this  section  we  summari/.c  the  relationships  wc  have  found  between  liglilning  flash  rates  and 
certain  weather  parameters.  The  details  of  the  derivations  can  be  found  in  Warber  and  Sinclair 
1 1995], 


To  be  useful,  these  weatlier  parameters  must  be  available  globally  and  forecast  up  to  several 
days  into  the  future;  or  else  tliey  must  be  derivable  from  sueh  parameters.  For  the  eoirclation 
study,  two  primary  data  sets  were  used;  a  lightning  flash  rate  archive  and  a  cloud 
piieiioiT'Ciiology  archive.  Archived  lightning  (lash  rate  tlata  is  available  from  the  National 
Lightning  Detection  Network  (NLDN),  a  system  of  100  sciisor.s  that  record  the  positittn,  time, 
number  of  return  suokes,  ;md  other  parameters  for  cloijd-li)-g,n)UiKi  fla.shes  in  the  U.S..  lire 
.sensor  sites  r.re  shown  in  iugurc  12-1.  Tlic  Nl.DN  claui  has  a  nominal  spatial  resolution  of  10 
kin  [Orville,  1990[.  Figure  12-2  shows  lightning  Hashes  deleeted  by  the  NLDN  in  four  UT  time 
blocks,  Cloiul  top  height  and  thickness  ((rplical  depth)  and  two  olhoi  parameters  were  derived 
fium  the  international  Satellite  Cloud  Climatology  IVoject  (ISCCP)  archive.  Ilie  ISCCFdata 
ineluiks  over  132  parameters  grouped  by  3-hour  time  blocks  and  a  nominal  250  km  spatial 
resolution.  The  correlation  study  utilizcil  data  from  both  NLDN  and  ISCCP  for  the  period  July 
1 990.  i'n  addition  to  the  ISCCP  data,  a  .set  of  Natitriial  Weather  Serviee  radar  summary  maps  for 


July  1990  and  a  set  of  weather  foreeast  maps  for  that  period  were  used  [Warber  and  Sinclair, 


1995[. 


Figure  12-1 .  Location  of  over  100  inagnetic  direcfloii  finders  in  the  National  Lightning  Detection 
Network  (NLDN). 
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Correlations  were  sought  between  lightning  flusli  rates  and  die  following  four  quantities 
derivable  from  tlie  ISCCP  observables: 


1.  Convective  cloud  top  height, 

2.  All  cloud  top  height, 

3.  Li*'' 'd  index,  L] 

4.  Convective  I  y  Available  Potential  Cineigy  (CAPC). 


In  addition,  correlation  between  tlasli  rates  and  a  heuristically  weighted  thunderstonn  density 
derived  by  analysis  of  weather  maps  was  studied.  'I’he  four  parameters  above  are  all  derivable 
by  standard  equations  from  the  basic  observables  recorded  in  the  ISCCP,  In  particular, 
convective  cloud  top  height  cati  be  inferred  Iry  discriminating  deep  convective  clouds  from 
others  using  the  ISCCP  cloud  optical  depth  ilata  in  combination  with  ISCCP  cloud  top  pre.ssure 
data.  Cloud  top  heights  are  computed  from  clouil  top  pressures  by  the  standard  calculation  for 
tire  altitude  at  which  an  ideal  gas  in  a  uniform  gravitational  field  reaches  a  given  prc:-;suie,  under 
the  a.s.suinption  that  the  rate  of  temperature  decrease  between  pre.ssure  levels  is  constant,  dire 
ISCe  data  set  also  contains  the  information  needed  to  calculate  the  lifted  index  {L[)  and  CAPE, 
both  of  which  are  related  to  tlie  virtual  temperatare  versus  allituile  profile  of  an  air  mass  (see 
Warber  and  Sinclair  [1995]  for  detailed  derivations). 


•  • 
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The  correlation  analysis  pruoeeiis  by  first  determining  a  least-squares  fit  between  flash  rate  and 
the  chosen  climatological  parameter  ifom  1  through  4  above.  Ali  the  relationships  could  be  fit 
as  least-squai’es  lines  on  log-log  or  log-linear  graphs.  The  fit  for  CAPlr  is  shown  in  Figure  12-3. 
Data  iTom  the  first  1 5  days  of  July  FJyO  were  used  to  deiermiiie  the  fits.  Data  iroints  with 
observed  flash  rates  of  zero  were  not  used.  After  a  fit  was  found,  data  from  the  second  15  days 
of  the  nu)nth  were  used  to  check  the  stability  of  the  fit.  During  tlie  check  all  data  points,  both 
zero  and  non -zero  ob.served  flash  rates  were  u.sed.  This  was  to  cheek  whether  tlic  fit,  which  is 
determined  from  only  non-zerti  flash  rates,  predicts  lightning  flashes  when,  in  fact,  there  arc 
none.  Two  correlation  coefficients  were  then  conijuited:  that  between  the  first  1 5  days  of  flash 
rate  data  and  the  fit  (.based  on  that  data)  and  that  between  the  second  1.5  days  of  tiata  and  the 
predicted  values  (using  the  fit  derived  from  the  first  15  days).  If  both  correlation  coefficients  are 
greater  than  0.5  (in  magnitude)  we  conclude  that  the  fit  is  a  stable  and  consistent  predictor  of 
flash  rate  from  the  given  weather  parameter.  In  instances  where  the  check  correlation  falls 
below  0.5,  however,  we  conclude  that  either  1)  the  slope  of  the  relation  between  flash  rate  and 
tile  cho.sen  weatlier  parameter  changes  significantly  over  the  course  of  a  mouth,  or  2)  tlie 
predictor  gives  significant  llitsh  rates  wlierc,  in  fact,  llicre  aie  none. 

In  particular,  the  lifted  indc,\  and  all  cloud  tti))  height,  Lj  and  were  fountl  to  he  inadequate 
predictors  of  llasli  rate  when  taken  by  themselves  (in  the  absence  of  otiier  ini'onnalion).  The 
following  conditions,  when  satisfied,  combined  L|  or  Z^., -based  predictions  raise  liic  check 
correlations  into  tlie  significance  range  (0.65  to  0.75): 

1)  If  there  is  a  thunderstorm  predicted  for  tiny  2.5"  x  2.5"  sub.square. 

2)  If  the  subsquarc  bad  lightning  in  it  during  the  previous  24  hours. 

3)  If  tJie  subsquare  is  over  a  mountainous  region  and  cloud  top  prc-ssure  is  less  Hum 
300  mb  (use  with  Lj  predictions) 

4)  If  the  subsquare  is  in  tlie  path  of  an  ongoing  .storm  system  (requires  radar  storm  summary 
maps). 

llicse  four  au.\iliaiy  conditions,  when  used  to  .select  r"..s-tri<'ted  subsets  of  weather  data  for 
lightning  flash  ate  predictions  will  be  referred  to  as  pmiictors  or  auxiliary  predictors. 

The  fits  between  flash  rate  and  the  four  weather  parameters  above,  Z^.,  Lj,  and  CAFF 


were: 


1)  Convective  cloud  top  height:  F  =  2.93  x  10'^'  .sec  *  .sq-*  (1  sq  s  5"  x  5“  area) 
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2)  All  cloud  top  height  (must  use  one  auxiliary  predictor):  F  =  3,73  x  10‘^  see*  sq'* 
i)  Lifted  Index  (must  use  one  auxiliary  predictor): 

P'nm  “  4.8  X  1()'2  exp(-0.42  Lj)  see*  sq"’  (non-niountain  region) 

F[^  =  2.7  X  lO"^- exp(-().47  Li)  sec"*  sq"*  (mountain  region) 

4)  CAPE: 

^NM  “  CAPE ’'25  .see*  sq-*  (non-mountain  regions) 

F[^  =  0.33  CAPE*  **^  .see*  .sq-*  (mountain  regions). 

A  fifth  prediction  metliod  derives  flash  rates  directly  from  the  extent  of  area  covered  by  severe 
weather.  A  few  days  worth  of  weatlier  maps  from  July  I  WO  were  digitized  and  overlaid  witii  a 
2.5“  X  2.5°  grid.  Each  subsiiuarc  containing  a  thunderstorm  prediction  was  marked  and  assigned 
tJie  value  1.  Then  each  marked  subsquarc  tliat  hud  marked  snbsquares  on  all  four  siilcs  was 
assigned  the  value  2.  Then  each  subsquare  with  a  value  2  that  had  adjacent  subsquares  with  the 
value  2  was  given  the  value  3.  The  process  was  repealeil  until  no  subsquares  were  changed. 

ITien  the  subsquares  were  assigned  the  following  flash  rates: 

•  Pnm  =  *’•2  X  value  *•’  •  n(h)  see"*  .sq'*  (non-mountain) 

•  Em  ~  value  ’-2  •  n(h)  .sec*  sq'*  (mountain) 

where  the  values  1.1  and  1.2  were  found  to  give  the  best  fit  to  observed  flash  rates  and  n(h)  is 
tire  normalized  local  time  dependence  factor  for  tlie  local  time  of  day  or  the  prediction.  Figure 
12  4  shows  the  factor  n(h).  (This  diurnal  curve  was  derived  by  grouping  all  the  July  1990 
NLDN  flash  rate  data  in  local  time  blocks.)  lire  overall  correlation  between  observation  and 
prediction  was  lower  than  that  using  tlie  other  methods.  In  the  ab.sence  of  otJier  data,  however, 
tlicse  weather  map  fits  can  be  used  to  make  rough  predictions. 

In  tire  following  we  di.seuss  tlie  lightning  forecasting  mcthiKls  we  have  developed  with  a  view 
toward  applications.  We  have  divided  the  lightning  forecasts  into  tliree  categories:  short,  mid, 
and  long  term  forecasts.  Short-term  forecasts  are  those  made  up  to  12  to  18  hours  in  advance, 
mid-term  forecasts  are  made  24  to  48  hours  in  advance,  aiul  long-term  forecasts  are  made  for 
periods  beyond  48  hours.  Tliese  types  of  forecasts  are  distinguished  by  the  methods  used  to 
make  them  and  the  confidence  wc  have  in  the  results. 


Figure  12-4.  Number  of  lighting  flashes  per  hour  as  a  function  of  local  time  in  July,  1990  over  all  NLDN  • 

sensors.  Normalized  to  1  at  1530  LT. 


12,2  Short-Term  Forecasts 

We  use  data  sources  in  order  of  reliability:  a)  direct  lightning  measurement  systems,  e.g., 

NLDN,  b)  cloud  top  heights  directly  from  weather  radar,  and  c)  cloud  top  heights  inferred  from 
satellite  images.  Using  a,  b,  or  c  we  detennine  the  current  lightning  distribution  (by  the  fits  for 
cloud  top  height  and  an  auxiliary  predictor  for  b,  c).  The  flash  rate  is  determined  for  each  2.5°  x 
2.5°  subsquare.  This  initial  distribution  is  called  the  flash  rate  basis.  Subsquares  in  the  basis  are 
then  grouped  into  storms.  The  storm  track  is  then  projected  up  to  24  hours  ahead  (using  the 
current  lighming  distributions  and  distributions  from  the  past  to  determine  how  the  center  of 
each  stomi  is  moving).  The  basis  flash  rates  are  copied  to  the  subsquares  around  the  new  storm 
center.  The  new  flash  rate  basis  is  adjusted  using  the  diurnal  curve  in  Figure  12-4.  If  we  write 
the  basis  flash  rate  as  and  the  br.sis  local  time  at  a  subsquare  as  to  then  the  adjusted  value  is: 

Fp  =  pQ  •  n(tp)/n(to) 


•  # 


» 


where  the  local  time  for  tlie  forecast  is  tp. 
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12.3  Mid-Term  Forecasts 


Mid-term  forecasts  will  be  based  either  on  predicted  CAPE  or  predicted  Lj  and  one  auxiliary 
predictor.  The  NWS  predicts  lifted  index  every  12  hours  for  he  periods  12,  36,  and  48  hours 
into  the  future. 

^2.4  Long-Term  Forecasts 

Long-term  forecasts  will  be  based  on  the  weather  map  correlation  (fifth  prediction  procedure 
abovef  This  requires  that  48  hour  severe  weather  predictions  are  available  regularly.  For  noise 
forecasts  based  on  long  term  forecasted  Hash  rates  it  is  more  reliable  to  give  ranges,  e.g.,  2  dB 
above  or  below  the  normal  average,  2  to  5  dB,  5  to  8  dB,  8  to  1 1  dB,  1 1  to  15  dB,  or  15  dB 
above  or  below  the  normal  average  noise  level  (determined  from  the  original  GLO  maps), 

12.5  Conclusions  and  Recommendations 

In  this  secdon  we  have  described  the  techniques  needed  to  create  an  operational  dynamic  noise 
prediction  system.  The  NoiseProp  noise  model  will  accept  DGLO  data  produced  by  the 
prediction  methods  discussed  above  to  generate  DGLO  files  which  arc  transparent  to  the 
NoiseProp  computer  code.  The  DGLO  prediction  codes  will  need  regular  access  to  data  from 
external  sources.  The  Navy's  Fleet  Numerical  Meteorology  and  Oceanographic  Center 
(FNMOC)  and  the  AFGWC  generate  the  kind  of  data  we  need  for  flash  rate  predictions  on  a 
global  scale.  AFGWC  has  expressed  interest  in  providing  forecast  data  which  would  be  needed 
in  the  lightning  prediction  codes.  It  seems  that  the  oest  way  to  implement  such  a  system  will  be 
for  a  central  location  to  collect  tire  needed  weather  data  and  process  it  into  the  DGLO  maps. 
These  maps  would  then  be  made  available  for  downloading.  A  computer  system  will  have  to  be 
created  to  completely  automate  the  forecasting  method  which  is  presently  far  too  labor  intensive 
to  operate  on  i  daily  basis.  The  forecast  algorithms  must  also  be  updated  over  time  as  larger 
numbers  of  predictions  and  measurements  are  correlated.  In  this  last  regard,  tlie  forecasting 
system  would  have  to  undergo  a  training  period.  The  present  report's  prediction  fits  are 
provisional  insofar  as  they  are  derived  fiom  one  month  of  daui  on  one  continent  only. 


13.0  REFERENCES 


Abramowitz,  M.  and  I.  A,  Stcgun,  Handbook  of  Mathematical  Functions,  National  ^Uireau  of 
Standards,  US  Government  Printing  Office,  Washington,  DC,  1972. 

Barghausen,  A.  F.  et  al.,  Predicting  Fong  Term  Operational  Parameters  of  High-Frequency 
Sky-Wave  Telecommunications  Systems,  ESSA  Technical  Report  ERL  1  lO-ITS  7S,  U.S. 
Department  of  Commerce,  Boulder,  CO,  1969. 

Bilitza,  D.  et  al..  International  Reference  Ionosphere  1990,  NSSDC/WDC-A-R&S  90-22, 
National  Space  Science  Data  CenterAVorld  Data  Center  A  for  Rockets  and  Satellites, 
Greenbelt,  MD,  November  1990. 

Bremmer,  H.  Terrestrial  Radio  Waves,  Eksevier,  Netherlands,  1949. 

Budden,  K.  B.,  Radiowaves  in  the  Ionosphere,  Cambridge  University  Press,  Cambridge,  UK, 
1966. 

- ,  The  Wave  Guide  Mode  Theory  of  Wave  Propagation,  LOGOS  Press,  London,  1961. 

Cambridge  Computation  Laboratory,  Tables  of  the  Modified  Hankel  Functions  of  Order  One- 
Third  and  of  Their  Derivatives,  Harvard  University  Pre.ss,  Cambridge,  MA,  1945. 

CCIR  (International  Radio  Consultative  Committee),  Report  340-1,  International 
Telecommunications  Union,  Geneva,  Switzerland,  1966. 

- ,  Report  340-2,  International  Telecommunications  Union,  Geneva,  Switzerland,  1967. 

. .  Interim  Method  for  Estimating  Sky  Wave  Field  Strength  and  Transmission  Loss  at 

Frequencies  Between  2  and  30  MHz,  Report  252-2,  International  Telecommunications  Union, 
Geneva,  Switzerland,  1986. 

. ,  World  Distribution  and  Characteristics  of  Atmospheric  Radio  Noise,  Report  322-3, 

International  Telecommunications  Union,  Geneva,  Switzerland,  1988. 

. ,  Sky-Wave  Propagation  and  Circuit  Perforrruince  at  Frequencies  Between  About  30  kHz 

arul  500  kHz,  Report  265-7,  International  Telecommunications  Union,  Geneva,  Switzerland, 
1990. 


— ,  CCIF^  Reference  Ionospheric  Characteristics  and  Methods  of  Basic  MIJF,  Operational 
MUr,  and  Ray-Path  Prediction,  Recommendation  434-.5,  International  Telecommunications 
Union,  Geneva,  Switzerland,  1992a. 


— ,  Sky-Wave  Field  Strength  Prediction  Method  for  the  Broadcasting  Service  in  the 
Frequency  Range  150  to  1600  kHz,  Recoimr-endation  435-7,  International 
Telecommunications  Union,  Geneva,  Switzerland,  iy92b. 


— ,  CCIR  HF  Propagation  Method,  Recommendation  533-3,  International  Telecommunica- 
tio.is  Union,  Geneva,  Switzerland,  1992c. 


Central  Radio  Propagation  Laboratory  (CRPL),  Ionospheric  Radio  Propagation,  National 
Bureau  of  Standards  Circular  462,  U.S.  Department  of  Commerce,  Washington,  DC,  June 
1948. 

Davies,  K.,  Ionospheric  Radio,  Peter  Peregrinus  Ltd,,  London,  United  Kingdom,  1990. 

Davis,  P.  M.,  Jr.,  Short-Term  Prediction  ofFl-Layer  Maximum  Usable  Frequencies  from  Local 
Magnetic  Activity,  Technical  Report  ERL  157-175-100,  Institute  for  Telecommunication 
Sciences,  Boulder,  CO,  1970. 

ITicker,  R.,  "A  Microcomputer  program  for  HF  Field  Strength  Prediction,"  Fifth  International 
Conference  on  Antennas  and  Propagation,  Part  2,  Propagation,  Institute  of  Electrical 
Engineers,  lEE,  United  Kingdom,  1987. 

Greifinger,  C.,  et  al..  Modeling  ofTM  and  TF,  VLF  Atmospheric  Noise,  Research  and 
Development  Associates.  RDA-TR-()226128902-001,  Augu.st  1989. 

Hatfield,  V  .  E.  and  G.  Smith,  AMBCOM  User's  Guide  for  Engineers,  Air  Force  Contract  FO 
8606-85-C-0018,  SRI  International,  Menlo  Park,  CA,  January  1987. 

Horner,  F.,  "Radio  Noise  From  niunderstorms,"  in  J.A.  Saxton  (Ed.j,  Advances  in  Radio 
Research,  Vol.  2,  Academic  Press,  New  York,  1964,  pp.  121-204, 

Horteiibach,  K.  J.,  and  F.  Rogler,  "On  the  Propagation  of  Short  Waves  Over  Very  Lung 
Distances:  Predictions  and  Observations,"  Telecommunicatiorus  Journal,  46,  pp.  320-327, 
1979. 

Inuki,  H.,  et  al.  "Development  of  the  Computer  Program  for  Predicting  LF  and  MF  .  Iky-Wave 
Field  Strengths,"  Radio  Research  Journal,  Vol.  29,  No.  152,  September  1983,  pp.  467-485, 
(in  Japane.se) 

Johler,  J.  R.,  "Propagation  of  the  Low-Frequency  Radio  Signal,"  Proceedings  of  the  IRE,  April, 
1962,  pp.  404-427. 


•  • 


-m-ti 


— ,  "On  the  Analysis  of  LF  Ionospheric  Radio  Propagation  Phenomena,"  Journal  of  Research 
of  the  National  Bureau  of  Standards -D.  Radio  Propagation,  Vol.  65D,  No.  5,  September' 
October  1961. 


■*f) 

(♦) 


Kotaki,  M.,  "Global  distribution  of  atmospheric  radio  noise  derived  from  Uiunderstorin 
activity,"  Journal  of  Atmospheric  andTerrestrial  Physics,  Vol.  46,  No.  10,  l‘^84,  pp.  867- 
877. 

Kotaki,  M.,  and  C.  Katoh,  "Global  Distribution  of  Atmospheric  Radio  Noi,se  Derived  from 
Distribution  of  Lightning  Activity,"  Journal  of  Radio  Research  Laboratory  (Tokyo),  Vol.  30, 
No.  129,  1983, 

LcVine,  D.  M.,  "Review  of  Measurements  of  the  RF  Spectrum  of  Radiation  from  Lightning," 
Meteorology  and  Atmospheric  Physics,  Vol.  37,  1987,  pp.  195-204. 

Lloyd,  J.  S.  e,t  al.,  Estimating  the  Performance  of  Telecommunications  Systems  Using  the 
Ionospheric  Transmission  Channel,  National  Telecommunications  and  Information 
Administration  (NTIA),  Boulder,  CO,  1978. 

Maxwell,  E.  L-.  et  al.,  Development  (fa  VLF  Atmospheric  Noise  Prediction  Model, 
Westinghouse  Geore.search  Laboratory,  Boulder,  CO,  Report  70- 1H2- VLF  NO-Rl,  30  June 
1970. 

Nanevicz,  J.  E.  et  al.,  "Observation  of  Lightning  in  the  Frequency  and  Time  Domains," 
Electromagnetics,  Vol.  7,  1987,  pp,  267-286. 

Orville,  R.  E.,  "Lightning  Ground  Flash  Density  in  the  Contiguous  United  States  -  1989," 
Monthly  Weather  Review,  Vol,  1 19,  No.  2.  1991,  pp.  573-577. 

PoKempner,  M.,  Comparison  of  Available  Methods  for  Predicting  Medium  Frequency  Sky- 
Wave  Field  Strengths,  NTIS  Report  80-42,  National  Telecommunications  and  Infomiation 
Administration,  Boulder,  CO,  1980. 

Sailors,  D.  B.  and  R.  B.  Rose,  FIF  Sky-Wave  Field  Strength  Predictions,  Technical  Report 
1624,  Naval  Command,  Control  and  Ocean  Surveillance  Center,  RDT&E  Division,  San 
Diego,  CA,  September  1993. 

Sailors,  D.,  A  Discrepancy  in  CCIR  Report  S22-3,  Radio  Noise  Model:  The  Probable  Cause 
and  a  Recommended  Solution,  I'echnical  Report  2496,  Naval  Command,  Control  and  Ocean 
Surveillance  Center,  RDT&E  Division,  San  Diego,  CA,  1993. 


i«) 

k 


•  9 


52 


X; 


Spaulding,  A.  D.  and  J.  S.  Washbuni,  Atmospheric  Radio  Noise:  Worldwide  Levels  and  Other 
Characteristics,  NTIA  Report  85-173,  US  Department  of  Coiirmerce,  Washington,  DC.  1985. 

Stuztman,  W.  L.  and  G.  A.  Thiele,  Antenna  Theory  and  Design,  John  Wiley,  New  York,  1987. 

Uman,  M.  A.,  The  Lightning  Discharge,  Academic  Press,  Inc.,  Orlando,  FL,  1987. 

Wait,  J.  R.  and  A.  M,  Conda,  "Pattern  of  an  Antenna  on  a  Curved  Lossy  Surface",  IRE 
Transactions  on  Antennas  and  Propagation,  AP-6,  1958,  pp.  348-360. 

Wait,  J.  R.,  "Diffractive  Corrections  to  the  Geometrii;al  Optics  of  Low  Frequency  Propagation," 
in  Desirant,  M.  and  J.  L.  Michiels  (Eds.)  Electromagnetic  Wave  Propagation,  Academic 
Press,  New  York,  1960,  pp.  87-101. 

- ,  Electromagnetic  Waves  in  Stratified  Media,  Pergamon  Press,  New  York,  1970. 

. ,  "Recent  Analytical  Investigations  of  Electromagnetic  Ground  Wave  Propagation  Over 

Inhomogeneous  Earth  Models,"  Proceedings  of  the  IEEE,  Vol.  62,  No.  8,  August  1974. 

Warber,  C.  R.  and  E.  C.  Field,  Jr.,  "A  Long  Wave  Transverse  Electric -Tran.sverse  Magnetic 
Noise  Prediction  Model,"  Radio  Science,  Vol.  30,  No.  3,  May-June  1995,  pp.  783-797. 

. ,  Long  Wave  Noise  Prediction,  Vol.  1 ,  Physical  Basis  of  the  Model,  Technical  Repoit 

DNA-TR-90-224-V1,  Defense  Nuclear  Agency,  Alexandria,  VA,  September  1991. 

Warber,  C.  R.  and  R.  F.  Sinclair,  Feasibility  of  Dynamic  Ixmgwave  Noise  Forecasting,  Report 
2605,  Pacific-Sierra  Research  Corporation,  Santa  Monica,  CA,  October  1995. 

Weidman,  C.  D.  and  E.  P.  Krider,  The  Amplitude  Spectra  of  Lightning  Radiation  Fields  in  the 
Interval  from  1  to  20  MHz,"  Radio  Science,  Vol.  21,  No.  6,  November-December  1986,  pp. 
964-970. 

- ,  "Lightning  Amplitude  Spectra  in  tlie  Interval  from  lOU  kHz  to  20  MHz,"  Geophysical 

Research  Letters,  Vol.  8,  No.  8,  August  1981,  pp.  931-934. 

Willet,  J.  C.,  et  al.,  "Lightning  Electromagnetic  Radiation  Field  Spectra  in  the  Interval  from  0.2 
to  20  MHz,"  Journal  of  Geophysical  Research,  Vol.  95,  No.  D12, ,  November  1990,  pp. 
20,367-20,387. 

Wrenn,  G.  L.  and  A.  3.  Rodger,  "Geomagnetic  Modification  of  the  Mid-Latitude  Ionosphere: 
Toward  a  Strategy  for  the  Improved  Forecasting  of  foF2,"  Radio  Science,  Vol.  24,  No.  1, 
January-February  1989,  pp.  99-11 1. 


53 


•  •  •  m  9 


o 


o 


Appendix 

THE  CARE  AND  FEEDING  OF  NOISEPROP 


1.  System  requirements 

NoiseProp  requires  Windows  95,  and  approximately  5  MB  of  disk  space. 

2.  Installation 

To  install  NoiseProp,  insert  disk  #1  in  the  drive.  From  the  Start  button,  select 
Run  a:\setup.  The  installation  program  will  create  a  Start  Menu  group  and 
program  icon. 


3.  Running  the  program. 

From  the  Start  button,  select  Programs  I  NoiseProp  I  NoiseProp. 

For  each  run,  a  time,  month,  frequency,  location,  bandwidth,  sunspot  number, 
and  exceedence  are  required.  An  optional  K-index  can  be  entered,  and  multiple 
frequencies  and  locations  can  entered. 

Most  parameters  are  entered  from  the  main  window.  Frequencies  and  locations 
are  entered  from  their  own  windows,  accessed  through  the  Edit  buttons  in  the 
main  window.  To  add  a  frequency  or  location,  enter  the  desired  value  and  click 
the  Add  button.  To  delete  a  value,  double-click  the  value  in  the  list-box,  and 
click  the  Delete  button.  To  modify  a  value,  double-click  the  value  in  the  list-box, 
enter  a  new  value  and  click  the  Replace  button.  To  finish  selecting  frequencies 
and/or  locations,  click  the  OK  button.  NOTE:  You  must  click  the  Add  button  to 
add  values  to  the  list.  Just  hitting  the  Enter  key  will  close  the  window  without 
adding  the  value  to  the  list. 

After  selecting  the  desired  parameters  for  a  run,  click  the  Run  button.  The 
generated  results  are  contained  in  the  file  RESULTS.OUT.  The  order  of  values 
in  the  file  is:  Frequency,  Time,  Longitude,  Latitude,  F,(Antenna  Noise  Figure),  F, 
(CCIR  prediction).  Units  for  F,  are  dB  above  kT^b,  where  k  is  Boltzmann’s 
constant,  T^  is  reference  temperature  (288  K),  and  b  is  receiver  bandwidth  (Hz). 

To  save  or  restore  parameters,  click  the  Save  or  Load  button,  and  enter  the 
desired  file  name. 


5.5 


